In [ ]:
More optimization trials¶
Previously, I tried throwing a bunch of optimizers at my model. I liked the results from the scipy.optimize routines COBYLA and BFGS. Let's try adjusting how we handle the constraints, and see if we get better results.
In [10]:
%load_ext autoreload
%autoreload 2
In [11]:
import matplotlib
import numpy
from IPython.display import HTML, SVG
In [12]:
#matplotlib.use('svg')
import matplotlib.pyplot as plt
plt.rcParams['svg.fonttype'] = 'none'
plt.rcParams['font.sans-serif'] = 'Arial'
In [13]:
from io import BytesIO
def pltsvg():
imgdata = BytesIO()
plt.savefig(imgdata)
imgdata.seek(0)
display(SVG(data=imgdata.read()))
In [14]:
import ammonia1
import system_aqua1
import scipy.optimize
In [19]:
T_heat_reject = 305.
UAgoal = 100
xB = [400, 1, T_heat_reject, 3, T_heat_reject, 5, 285, 4, T_heat_reject, 0.15]
bdry = system_aqua1.makeBoundary(xB)
Trial 2.1¶
In general, here is how we compute the modified objective function. GenOpt manual suggests adding functions to the objective, namely separate barrier (B) and penalty (P) functions. The objective function is, with step number as $k$,
$$ Objective = -Q_{cooling}(x) + \mu_1 B(x) + \mu_2 P(x) $$$$ \mu_1 = k^{-2}$$$$ B(x) = \left(\sum_{j=1}^{N_B} g_j(x) \right) ^ {-1} $$$$\mu_2 = k^{2}$$$$ P(x) = \sum_{j=N_B+1}^{N_B + N_P} g_j(x)$$GenOpt implements a mechanism to pass in the step number. However, using the scipy.optimize routines, passing in the step number is a challenge, so the Problem object should probably use its own counter of function calls. This could be quirky, since then the objective function is changing with each call. Some routines allow for termination after a number of steps, so we could run a while, increment step number, then continue.
But we need to make some decisions.
- Which constraints should be handled by barriers, and which by penalties?
- Use barriers for hard constraints, eg. required for feasibility of the chiller model
- Use penalties for soft constraints, eg. heat exchange feasibility that can be violated now and satisfied later
- With scipy.optimize, we can directly apply constraints on the input variables (such as $T_1 > T_2$)
In [15]:
def saturate(x, bottom=-numpy.inf, top=0):
a_bottom = numpy.empty_like(x)
a_top = numpy.empty_like(x)
a_bottom.fill(bottom)
a_top.fill(top)
return numpy.minimum(a_top,
numpy.maximum(a_bottom,
x))
def barrier1(c, length_scale):
"""The domain for B is the feasible set only.
We want B large near boundary, but small elsewhere.
Feasible means c > 0 and g < 0. Please input c.
"""
g = numpy.negative(c)
b = numpy.exp(g / length_scale)
return numpy.sum(b)
def decay1(step_number, initial_value = 1., rate = 1.):
"""A decaying function to scale barrier functions.
step_number: as this increases, I decay.
initial_value: value when step_number = 0.
decay_rate: how fast to decay.
"""
# Magnitude tends to zero, slowly
mu_B = initial_value * numpy.exp(-rate * step_number)
return mu_B
def penalty1(c, length_scale):
"""We want P = 0 for feasible, P > 0 elsewhere
Feasible means c > 0 and g < 0.
"""
g = numpy.negative(c)
g = saturate(g,bottom=0,top=numpy.inf)
p = (g / length_scale)**2
return numpy.sum(p)
def grow1(step_number, initial_value = 1., rate = 1.):
"""A growing function to scale penalty functions."""
# Magnitude tends to infinite, bit more quickly
mu_P = initial_value * numpy.exp(rate * step_number)
return mu_P
In [41]:
# Demo barriers
xx = numpy.linspace(0,1)
yy = numpy.linspace(0,1)
X, Y = numpy.meshgrid(xx,yy)
B = numpy.zeros_like(X)
mu_B = 10
for i, (xi,yi) in enumerate(zip(X.flat,Y.flat)):
# The domain for B is the feasible set only.
# We want B large near boundary, but small elsewhere.
# Feasible means c > 0 and g < 0.
c = numpy.array([xi, yi, (1-xi), (1-yi)])
# The length scale for the barrier
dx = 0.05
B.flat[i] = mu_B * barrier1(c,dx)
plt.close('all')
fig = plt.figure()
plt.plot(X[25],B[25],'r')
plt.title("barrier1()")
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X,Y,B,cmap=cm.coolwarm)
ax.plot(X[25],Y[25],B[25],'r')
#pltsvg()
plt.show()
xx = numpy.linspace(-1,11)
yy = numpy.linspace(-1,11)
X, Y = numpy.meshgrid(xx,yy)
P = numpy.zeros_like(X)
mu_P = 10
for i, (xi,yi) in enumerate(zip(X.flat,Y.flat)):
# We want P = 0 for feasible, P > 0 elsewhere
# Feasible means c > 0 and g < 0.
c = numpy.array([xi, yi, (10-xi), (10-yi)])
# The length scale for the penalty
dx = 1
P.flat[i] = mu_P * penalty1(c, dx)
plt.close('all')
fig = plt.figure()
plt.plot(X[25],P[25])
plt.title("penalty1")
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X,Y,P,cmap=cm.coolwarm)
ax.plot(X[25],Y[25],P[25],'b')
#pltsvg()
plt.show()
In [44]:
step_number = numpy.arange(50)
mu_B = decay1(step_number,initial_value=1000,rate=0.2)
mu_P = grow1(step_number,initial_value=1,rate=0.2)
plt.close('all')
plt.figure()
plt.semilogy(step_number,mu_B,'r',label='barrier')
plt.semilogy(step_number,mu_P,'b',label='penalty')
plt.xlabel('Optimization step number')
plt.ylabel('Scaling for constraint function')
plt.legend()
plt.show()
In [21]:
step_number = 39
mu_B = 1 * numpy.exp(0.1 * step_number)
mu_B
Out[21]:
In [181]:
class Problem_2_1:
def __init__(self, bdry, UAgoal, mu=0.1):
self.bdry = bdry
self.UAgoal = UAgoal
self.mu = mu
self.Ncons = 7
self.n_calls = 0
# Soft constraints mode: this is sent to minimizer
self.constraints = [{'type': 'ineq',
'fun': self.constraint,
'args': (i,)
} for i in range(self.Ncons)]
def objective(self, xC):
step_number = numpy.floor(self.n_calls / 7)
self.n_calls += 1
print(xC,flush=True)
Q,B,P = 0.,0.,0.
try:
ch = system_aqua1.makeChiller(xC)
sys = system_aqua1.System(self.bdry, ch)
# Barriers
# Magnitude tends to zero, slowly
mu_B = 1000. * numpy.exp(-0.1 * step_number)
length_scale_b = 1
# Or ... magnitude fixed, but shape changes
mu_B = 1000
length_scale_b = 1 * numpy.exp(-0.1 * step_number)
# These are zero at the boundary ...
barriers = [ch.check_rectifier_delta_T] \
+ [deltaT
for name, deltaT, epsilon, UA, Qhx in sys.data]
B = mu_B * barrier1(barriers,length_scale_b)
# Penalties
# Magnitude tends to infinite
mu_P = 1 * numpy.exp(0.1 * step_number)
penalties = [(1.0 - sys.totalUA/self.UAgoal)]
P = mu_P * penalty1(penalties,1)
Q = sys.chiller.Q_evap
except:
Q = numpy.inf
print(self.n_calls, step_number, Q, B, P, "\n", flush=True)
return -Q + B + P
def constraint(self, x, *args):
cons = [x[0] - 0.1,
1. - x[0],
x[2] - x[1] - 1.0,
x[3] - x[2] - 0.1,
x[4] - x[1] - 10.0,
x[5] - x[3] - 1.0,
x[5] - x[4] - 1.0]
if len(args) > 0:
i, = args
return cons[i]
else:
return cons
def callback(self, x):
print("Did an iteration at ", x)
In [ ]:
P = Problem_2_1(bdry, UAgoal)
In [183]:
rT = T_heat_reject
x = numpy.array([0.05, 278.45, rT+7, rT+8, rT+5, 395.15])
In [125]:
x = numpy.array([0.51284472, 277.97717012, 312.16427764, 313.6952877,
310.24856734, 374.14020482])
In [205]:
P.n_calls = 7 * numpy.floor(P.n_calls / 7)
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints, callback=P.callback,
options={'disp':True,'maxiter':1000,'rhobeg':0.05})
opt
Out[205]:
In [206]:
x = opt.x
x
Out[206]:
In [ ]:
# Let's bebug bad points encountered during optimization.
# For reference, here is a good case.
good_X1 = numpy.array([0.51284472,
277.97717012,
312.16427764,
313.6952877,
310.24856734,
374.14020482])
# Now, let's compare some degenerate cases.
# This one requires absorber stream to compute points at high concentrations.
bad_X1 = numpy.array([0.2,
290.95402556,
306.64229519,
313.6952877,
300.86211291,
384.36256216])
# This one is not feasible, insofar as it yields positive rectifier heat flow,
# that is, heat *input* to rectifier, which should be rejected as "bad" design.
# The cause is that T_abs is very low, so the rich stream is very rich,
# and in fact the rectifier isn't needed...
# So, I'm adding a constraint to keep the absorber warmer than evaporator.
bad_X2 = numpy.array([1.,
293.03582607,
308.98651075,
313.6952877,
293.29251019,
386.99474])
# This one is also degenerate in the rectifier.
# Try increasing the padding on the constraints.
bad_X3 = numpy.array([1.,
292.26696584,
301.95207453,
314.19203364,
293.26696584,
389.23622472])
bad_X4 = numpy.array([1.,
277.22895294,
278.22895294,
308.37781922,
282.22895294,
417.08853392])
x = good_X1
In [207]:
ch = system_aqua1.makeChiller(x)
display(ch)
P.constraint(x)
Out[207]:
In [208]:
sys = system_aqua1.System(bdry, ch)
sys
Out[208]:
In [209]:
1-sys.totalUA/P.UAgoal
Out[209]:
In [210]:
penalty1([(1 - sys.totalUA/P.UAgoal)], 1)
Out[210]:
In [211]:
box = numpy.ones_like(x)
box[0] = 0.05
ss = numpy.linspace(-1,1,5)
yyy = numpy.zeros([len(box),len(ss)])
for i,boxi in enumerate(box):
yy = numpy.zeros_like(ss)
for j, s in enumerate(ss):
x_vary = x.copy()
x_vary[i] += boxi * s
yyy[i,j] = P.objective(x_vary)
if (numpy.array(P.constraint(x_vary)) < 0).any():
yyy[i,j] = numpy.nan
In [213]:
plt.close('all')
plt.figure(figsize=(10, 4))
ax1=plt.subplot(1,len(box), 1)
for i,_ in enumerate(box):
ax=plt.subplot(1, len(box), i+1, sharey=ax1)
xxx = x[i] + ss * box[i]
plt.plot(xxx, yyy[i,:], '.-')
plt.setp(ax.get_yticklabels(), visible=bool(i==0))
plt.xlim(min(xxx),max(xxx))
plt.ylim(min(yyy.flat),min(yyy.flat)*0.9)
plt.show()
In [173]:
x = numpy.array([ 0.31, 277.8, 313.,
314.2, 311.5, 394.5])
Trial 2.2¶
Just brainstorming ... As a brute-force backup, I would like to have an optimize routine that is a gradient-based method where the gradient calculation and line-search algorithms simply react to passing over the feasible boundary using exceptions instead of constraint functions. (This matters because the total UA constraint function can only be computed when the feasible heat exchange constraints are already satisfied, so in some cases it would have to return infinite values.) In other words, this would be an adaptive step size method. In pseudo-code,
In [ ]:
def optimize(fun,x0,max_iter = 100):
x = x0.copy()
for iter in range(max_iter):
if step == 0:
break
G, H = calcGradient(fun, x)
dx = calcLineSearch(fun, x, G, H)
x = x + dx
f = fun(x)
return x, f
def calcGradient(fun, x0, max_iter = 100):
max_iter = 100
# Initial step sizes by component
deltax = ones_like(x0)
# Function values to the right
fplus = zeros_like(x0)
# Function values to the left
fminus = zeros_like(x0)
f0 = fun(x0)
for j, _ in enumerate(deltax):
for i in range(max_iter):
dx = numpy.diag(deltax)[j]
try:
xplus = x0 + dx
xminus = x0 - dx
fplus[j] = fun(xplus)
fminus[j] = fun(xminus)
break
except:
deltax[j] *= 0.5
grad, hess = calcGradHelper(f0, deltax, fplus, fminus) # TODO
return grad, hess
def calcLineSearch(fun, x0, G, H):
max_iter = 100
step = 0
direction = lshelper1(x0, G, H)
for i in range(max_iter):
try:
lshelper(f
except:
step *= 0.5
return step * direction
Aside¶
Furthermore, here is another related strategy. I could optimize first without the constraint on total UA. That would give me a feasible system. Then, given the desired total UA, I could have a feasible starting point by then turning down the mass flow rate using a linesearch on that parameter.
Trial 2.3¶
Another idea (suggested in documention of early trials) is to do something different with the heat exchange feasibility constraints, utilizing energy imbalance. The final goal would be to allow UA values as inputs, and optimize the cooling capacity, as intended. An intermediate goal may be to simply demonstrate a solver with UA values as inputs (using some optimization routine to determine the temperature points that achieve the specified UA values).
Trial 2.4¶
New idea: use machine learning to build constraint functions, as follows. Let the optimization routine run. For each input evaluated, classify its feasibility (bonus: add dimensions for root cause of failure). Gather the data for classified points and generate a classification function through regression, etc. Restart the optimizer with this function
Trial 2.5¶
An extension of a previous idea, kind of the dual problem -- Perform optimization in two steps: first with no constraint on UA (cost) to determine maximum Q (output). Then for each Q in a set of lesser values, optimize with respect to cost. This will deliver an optimal cost vs output curve.
For the first step, since we are not interested in UA, we can be less strict with the constraints on heat exchange feasibility. Specifically, we can impose them via penalty constraints instead of barrier constraints.
In [92]:
class Problem_2_5_A:
def __init__(self, bdry, UAgoal, mu=0.1):
self.bdry = bdry
self.UAgoal = UAgoal
self.mu = mu
self.Ncons = 7
self.n_calls = 0
# Soft constraints mode: this is sent to minimizer
self.constraints = [{'type': 'ineq',
'fun': self.constraint,
'args': (i,)
} for i in range(self.Ncons)]
def objective(self, xC):
step_number = numpy.floor(self.n_calls / 7)
self.n_calls += 1
print(xC,flush=True)
Q,B,P = 0.,0.,0.
try:
ch = system_aqua1.makeChiller(xC)
sys = system_aqua1.System(self.bdry, ch)
# Barriers
# Magnitude tends to zero, slowly
mu_B = 1000. * numpy.exp(-0.1 * step_number)
length_scale_b = 1
# Or ... magnitude fixed, but shape changes
mu_B = 1000
length_scale_b = 1 * numpy.exp(-0.1 * step_number)
# These are zero at the boundary ...
barriers = [ch.check_rectifier_delta_T]
B = mu_B * barrier1(barriers,length_scale_b)
# Penalties
# Magnitude tends to infinite
mu_P = 1 * numpy.exp(0.3 * step_number)
penalties = [deltaT - 0.01
for name, deltaT, epsilon, UA, Qhx in sys.data]
P = mu_P * penalty1(penalties,1)
Q = sys.chiller.Q_evap
except:
Q = numpy.inf
print(self.n_calls, step_number, Q, B, P, "\n", flush=True)
return -Q + B + P
def constraint(self, x, *args):
cons = [x[0] - 0.1,
1. - x[0],
x[2] - x[1] - 1.0,
x[3] - x[2] - 0.1,
x[4] - x[1] - 10.0,
x[5] - x[3] - 1.0,
x[5] - x[4] - 1.0]
if len(args) > 0:
i, = args
return cons[i]
else:
return cons
def callback(self, x):
print("Did an iteration at ", x)
In [93]:
P = Problem_2_5_A(bdry, UAgoal)
rT = T_heat_reject
x = numpy.array([0.05, 278.45, rT+7, rT+8, rT+5, 395.15])
In [119]:
P.n_calls = 7 * numpy.floor(P.n_calls / 7)
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints, callback=P.callback,
options={'disp':True,'maxiter':1000,'rhobeg':0.001})
opt
Out[119]:
In [120]:
x = opt.x
x
Out[120]:
In [121]:
ch = system_aqua1.makeChiller(x)
display(ch)
In [122]:
P.constraint(x)
Out[122]:
In [123]:
sys = system_aqua1.System(bdry, ch)
sys
Out[123]:
In [124]:
box = numpy.ones_like(x)
box[0] = 0.05
ss = numpy.linspace(-1,1,5)
yyy = numpy.zeros([len(box),len(ss)])
for i,boxi in enumerate(box):
yy = numpy.zeros_like(ss)
for j, s in enumerate(ss):
x_vary = x.copy()
x_vary[i] += boxi * s
yyy[i,j] = P.objective(x_vary)
if (numpy.array(P.constraint(x_vary)) < 0).any():
yyy[i,j] = numpy.nan
In [125]:
plt.close('all')
plt.figure(figsize=(10, 4))
ax1=plt.subplot(1,len(box), 1)
for i,_ in enumerate(box):
ax=plt.subplot(1, len(box), i+1, sharey=ax1)
xxx = x[i] + ss * box[i]
plt.plot(xxx, yyy[i,:], '.-')
plt.setp(ax.get_yticklabels(), visible=bool(i==0))
plt.xlim(min(xxx),max(xxx))
plt.ylim(min(yyy.flat),min(yyy.flat)*0.9)
plt.show()
Now that we have a guess for the maximum cooling we might get, let's optimize the cost function instead, and ask for that much cooling using a constraint. I'm thinking that a penalty constraint might be good, so that we add this to the objective:
$$ P = \left(Q - Q_{goal}\right) ^ 2 $$In [134]:
Q_opt = ch.Q_evap
Q_opt
Out[134]:
In [140]:
Q_goal = 0.9 * Q_opt
Q_goal
Out[140]:
In [170]:
class Problem_2_5_B:
def __init__(self, bdry, Q_goal, mu=0.1):
self.bdry = bdry
self.Q_goal = Q_goal
self.mu = mu
self.Ncons = 7
self.n_calls = 0
# Soft constraints mode: this is sent to minimizer
self.constraints = [{'type': 'ineq',
'fun': self.constraint,
'args': (i,)
} for i in range(self.Ncons)]
def objective_raw(self, xC):
try:
ch = system_aqua1.makeChiller(xC)
sys = system_aqua1.System(self.bdry, ch)
UA = sys.totalUA
except:
UA = numpy.nan
return UA
def objective(self, xC):
step_number = numpy.floor(self.n_calls / 7)
self.n_calls += 1
print(xC,flush=True)
UA,B,P = 0.,0.,0.
try:
ch = system_aqua1.makeChiller(xC)
sys = system_aqua1.System(self.bdry, ch)
# Barriers
# Magnitude tends to zero, slowly
mu_B = 1000. * numpy.exp(-0.1 * step_number)
length_scale_b = 1
# Or ... magnitude fixed, but shape changes
mu_B = 1000
length_scale_b = 1 * numpy.exp(-0.1 * step_number)
# These are zero at the boundary ...
barriers = [ch.check_rectifier_delta_T] \
+ [deltaT
for name, deltaT, epsilon, UA, Qhx in sys.data]
B = mu_B * barrier1(barriers,length_scale_b)
# Penalties
# Magnitude tends to infinite
mu_P = 1 * numpy.exp(0.3 * step_number)
penalties = [ch.Q_evap - self.Q_goal]
P = mu_P * penalty1(penalties,1)
UA = sys.totalUA
except:
UA = numpy.inf
print(self.n_calls, step_number, UA, B, P, "\n", flush=True)
return UA + B + P
def constraint(self, x, *args):
cons = [x[0] - 0.1,
1. - x[0],
x[2] - x[1] - 1.0,
x[3] - x[2] - 0.1,
x[4] - x[1] - 10.0,
x[5] - x[3] - 1.0,
x[5] - x[4] - 1.0]
if len(args) > 0:
i, = args
return cons[i]
else:
return cons
def callback(self, x):
print("Did an iteration at ", x)
In [171]:
P = Problem_2_5_B(bdry, Q_goal)
rT = T_heat_reject
x2 = x.copy()
In [159]:
P.n_calls = 7 * numpy.floor(P.n_calls / 7)
opt = scipy.optimize.minimize(P.objective, x2, method="COBYLA",
constraints=P.constraints, callback=P.callback,
options={'disp':True,'maxiter':100,'rhobeg':0.01})
opt
Out[159]:
In [182]:
x2 = opt.x
x2
Out[182]:
In [161]:
ch2 = system_aqua1.makeChiller(x2)
display(ch2)
In [162]:
P.constraint(x2)
Out[162]:
In [163]:
sys2 = system_aqua1.System(bdry, ch2)
sys2
Out[163]:
In [183]:
box = numpy.ones_like(x2)
box[0] = 0.05
ss = numpy.linspace(-1,1,5)
yyy = numpy.zeros([len(box),len(ss)])
for i,boxi in enumerate(box):
yy = numpy.zeros_like(ss)
for j, s in enumerate(ss):
x_vary = x2.copy()
x_vary[i] += boxi * s
yyy[i,j] = P.objective_raw(x_vary)
if (numpy.array(P.constraint(x_vary)) < 0).any():
yyy[i,j] = numpy.nan
In [184]:
yyy
Out[184]:
In [187]:
plt.close('all')
plt.figure(figsize=(10, 4))
ax1=plt.subplot(1,len(box), 1)
for i,_ in enumerate(box):
ax=plt.subplot(1, len(box), i+1, sharey=ax1)
xxx = x2[i] + ss * box[i]
plt.plot(xxx, yyy[i,:], 'o-')
plt.setp(ax.get_yticklabels(), visible=bool(i==0))
plt.xlim(min(xxx),max(xxx))
#plt.ylim(min(yyy.flat)*0.9,min(yyy.flat))
plt.show()
The sensitivity plot could go further by overlaying the constraint functions. Additionally, it would be nice to find the eigenvectors of the objective function and trace in those directions, indicating how far it is to the next constraint. But that is icing on the cake, so meanwhile I will continue with this approach in a new notebook and generate some parametric study data.
Scratch work¶
In [15]:
a=5
for i in range(9999):
a*=0.5
if a==0:
break
i,a
Out[15]:
In [25]:
gradient_step = numpy.array([1.,2.,3.,4.])
for j, _ in enumerate(gradient_step):
for i in range(2):
gradient_step[j] *= 0.5
gradient_step
Out[25]:
In [30]:
numpy.diag(gradient_step)[3]
Out[30]:
In [7]:
%%html
<pre id="TOC-markdown">TOC will be here</pre>
<script>
$("#TOC-markdown").html(
$('h1,h2,h3,h4').filter(":has(a)").map(function(){return " ".repeat($(this).prop("tagName")[1]) + "- ["+$(this).text()+"](" + $(this).children().attr("href") + ")";}).get().join("\n")
);
</script>
In [ ]: