HGDL Constrained Optimization of Schwefel’s Function
In this script, we show how HGDL is used for constrained optimization. Unconstrained optimization is simpler and automatically included.
#!pip install hgdl==2.1.9
First some function to make nice plots
%load_ext autoreload
%autoreload 2
import numpy as np
import plotly.graph_objects as go
def plot(x,y,z,data = None, constr = None):
fig = go.Figure()
fig.add_trace(go.Surface(x = x, y = y,z=z))
if data is not None:
fig.add_trace(go.Scatter3d(x=data[:,0], y=data[:,1], z=data[:,2] + 50,
mode='markers'))
if constr is not None:
fig.add_trace(go.Scatter3d(x=constr[:,0], y=constr[:,1], z=constr[:,2],
mode='markers'))
fig.update_layout(title='Surface Plot', autosize=True,
width=800, height=800, font=dict(
family="Courier New, monospace",
size=18),
margin=dict(l=65, r=50, b=65, t=90))
fig.show()
def make_plot(data = None, constraint1 = None, constraint2 = None):
x1,x2 = np.linspace(-500,500,100),np.linspace(-500,500,100)
x_pred = np.transpose([np.tile(x1, len(x2)), np.repeat(x2, len(x1))])
r1 = np.sqrt(160000.)
r2 = np.sqrt(40000.)
c1,c2 = r1*np.cos(np.linspace(0,2.*np.pi,100)),r1*np.sin(np.linspace(0,2.*np.pi,100))
d1,d2 = r2*np.cos(np.linspace(0,2.*np.pi,100)),r2*np.sin(np.linspace(0,2.*np.pi,100))
c3 = np.zeros((len(c2)))
d3 = np.zeros((len(c2)))
x1,x2 = np.meshgrid(x1,x2)
z = np.zeros((10000))
zc1 = np.zeros((10000))
zc2 = np.zeros((10000))
for i in range(10000):
z[i] = schwefel(x_pred[i], 1, 1)
if constraint1: zc1[i] = constraint1(x_pred[i])
if constraint2: zc2[i] = constraint2(x_pred[i])
for i in range(len(c1)):
c3[i] = schwefel(np.array([c1[i],c2[i]]))
d3[i] = schwefel(np.array([d1[i],d2[i]]))
plot(x1,x2,z.reshape(100,100).T, data = data,
constr = np.row_stack([np.column_stack([c1,c2,c3]),np.column_stack([d1,d2,d3])]))
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[2], line 4
2 get_ipython().run_line_magic('autoreload', '2')
3 import numpy as np
----> 4 import plotly.graph_objects as go
5 def plot(x,y,z,data = None, constr = None):
6 fig = go.Figure()
ModuleNotFoundError: No module named 'plotly'
Defining the Constraints and some Bounds
Keep in mind that not all local optimizers allow any combination of bounds and constraints Visit https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html for more information on that
bounds = np.array([[-500,500],[-500,500]])
def g1(x):
return (np.linalg.norm(x)**2/10.0) - 16000.0
def g2(x):
return (np.linalg.norm(x)**2/10.0) - 4000.0
Now we import HGDL and run a constrained optimization
from hgdl.hgdl import HGDL as hgdl
from hgdl.support_functions import *
import time
#example arguments that will be passed to the function, the gradient and the Hessian function
a = 5
b = 6
#constraint definitions form scipy
from scipy.optimize import NonlinearConstraint
nlc = NonlinearConstraint(g1, -np.inf, 0)
nlc = NonlinearConstraint(g2, 0,np.inf)
a = hgdl(schwefel, schwefel_gradient, bounds,
hess = None, ##if this is None, the Hessian will be approximated if the local optimizer needs it
global_optimizer = "random", #there are a few options to choose from for the global optimizer
#global_optimizer = "genetic",
local_optimizer = "dNewton", #dNewton is an example and will be changed automatically to "SLSQP" because constraints are used
number_of_optima = 30000, #the number fo optima that will be stored and used for deflation
args = (a,b), num_epochs = 1000, #the number of total epochs. Since this is an asynchronous algorithms, this number can be very high
constraints = (nlc,) #the constraints
)
a.optimize(x0=None)
res = a.get_latest()
for entry in res: print(entry)
res = a.kill_client()
Making a Plot
You should see the constraints and the foudn optima. If everything worked, the found points are in between the two constraints.
data = [np.append(entry["x"],entry["f(x)"]) for entry in res]
make_plot(data = np.array(data), constraint1 = g1, constraint2 = g2)