Trial 2.3 - Use inputs that are more physical¶
The idea here is to do something different with the heat exchange feasibility constraints. In previous trials, we have written a solver that inputs temperatures to solve the chiller cycle; then heat exchange feasibility is achieved by handing various constraints to the optimizer.
Here we want to solve a system with UA inputs. One method to accomplish this can be easily implemented by using the existing code and a numerical solver as follows. With a guess for T inputs, we solve the system and for each heat exchanger, we can build the (T,q) curves on both hot and cold sides. Then we evaluate Q resulting from the given UA and compare it to Q determined by the chiller cycle. The vector of [disparity in Q for each heat exchanger] will be minimized by varying T. However, it should be noted that the numerical solver adds a computational burden, and furthermore my existing method to computate Q given UA also uses a numerical solve. So in short, although the choice of inputs may remove the need for some temperature constraints, it would yield a slow simulation.
In [1]:
%load_ext autoreload
%autoreload 2
In [2]:
import matplotlib
import numpy
from IPython.display import HTML, SVG
In [3]:
#matplotlib.use('svg')
import matplotlib.pyplot as plt
plt.rcParams['svg.fonttype'] = 'none'
plt.rcParams['font.sans-serif'] = 'Arial'
In [4]:
from io import BytesIO
def pltsvg():
imgdata = BytesIO()
plt.savefig(imgdata)
imgdata.seek(0)
display(SVG(data=imgdata.read()))
In [5]:
import ammonia1
import system_aqua1
import scipy.optimize
In [6]:
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)
In [7]:
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 [49]:
type(numpy.sum([1,2,3]))
Out[49]:
In [50]:
class Problem_2_3:
def __init__(self, bdry, UAs):
"""Note: UAs to input are like this:
[gen,rect,abs,cond,evap]
"""
self.bdry = bdry
self.UAs = UAs
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
cost,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.1 * step_number)
penalties = []
P = mu_P * penalty1(penalties,1)
Q = numpy.zeros_like(sys.Q)
Q_excess = numpy.zeros_like(sys.Q)
for i, (name, hx) in enumerate(sys.hxs.items()):
Q[i] = hx.calcQ(self.UAs[i])
Q_excess[i] = Q - sys.Q[name]
print(Q)
cost = numpy.sum(numpy.square(Q_excess))
#except Exception as e:
# print(e)
# cost = numpy.inf
print(self.n_calls, step_number, cost, B, P, "\n", flush=True)
return cost + 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 [51]:
# [gen,rect,abs,cond,evap]
P = Problem_2_3(bdry, [30,1,30,30,30])
In [52]:
rT = T_heat_reject
x = numpy.array([0.5, 278.45, rT+7, rT+8, rT+5, 395.15])
In [53]:
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':100,'rhobeg':0.1})
opt
In [33]:
x = opt.x
x
Out[33]:
In [34]:
ch = system_aqua1.makeChiller(x)
display(ch)
P.constraint(x)
Out[34]:
In [35]:
sys = system_aqua1.System(bdry, ch)
sys
Out[35]:
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 [ ]: