Optimization method 2.5 parametric studies¶
Previously, I obtained reasonable convergence on optima using a two-step strategy: first, simply enforce feasibility and optimize cooling output without regard to cost; then for any given cooling output goal less than the maximum, enforce that as a constraint and optimize system cost (from size and flow parameters). With some art in selecting constraint functions, I implemented this using the COBYLA optimization routine.
Since that method can be made to work, let's do a parametric study in two parameters: temperature for heat rejection (T) and cost (UA). For each heat reject temperature, we will determine the maximum cooling (Q_max), then build a cost vs cooling curve. The data should form a surface of points (T,UA,Q). We can also visualize that surface by projecting and plotting isolines, for example, isolines of UA on a (T,Q) plot.
Setup¶
In [14]:
%load_ext autoreload
%autoreload 2
In [15]:
import matplotlib
import numpy
from numpy import array, inf, nan
import pandas
from IPython.display import HTML, SVG, clear_output
In [16]:
#matplotlib.use('svg')
%matplotlib notebook
import matplotlib.pyplot as plt
plt.rcParams['svg.fonttype'] = 'none'
plt.rcParams['font.sans-serif'] = 'Arial'
In [17]:
from io import BytesIO
def pltsvg():
imgdata = BytesIO()
plt.savefig(imgdata)
imgdata.seek(0)
display(SVG(data=imgdata.read()))
In [18]:
import ammonia1
import system_aqua1
import scipy.optimize
In [19]:
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 [16]:
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()
Strategy¶
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 [8]:
class Problem_2_5_A:
def __init__(self, bdry, mu=0.1):
self.bdry = bdry
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 KeyboardInterrupt as e:
raise e
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 [9]:
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 KeyboardInterrupt as e:
raise e
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)
Parametric study across T_heat_reject¶
In [20]:
def display_result(x):
try:
print("X = ")
display(x)
ch = system_aqua1.makeChiller(x)
sys = system_aqua1.System(bdry, ch)
display(sys)
except KeyboardInterrupt as e:
raise e
except Exception as e:
display(e)
In [21]:
t_heat_reject_range = 273.15 + numpy.arange(20,61,5)
In [147]:
opts = [None] * len(t_heat_reject_range)
Step one¶
In [24]:
for i,T_heat_reject in enumerate(t_heat_reject_range):
if opts[i] is None:
print("Case {}, T_heat_reject = {} K".format(i, T_heat_reject))
rT = T_heat_reject
xB = [400, 1, T_heat_reject, 3, T_heat_reject, 5, 285, 4, T_heat_reject, 0.15]
bdry = system_aqua1.makeBoundary(xB)
P = Problem_2_5_A(bdry)
x = numpy.array([0.05, 278.45, rT+7, rT+8, rT+5, 395.15])
P.n_calls = 7*30
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints, callback=P.callback,
options={'disp':True,'maxiter':50,'rhobeg':0.1})
x = opt.x
display("First pass and we have this ... ")
display_result(x)
P.n_calls = 7*40
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints, callback=P.callback,
options={'disp':True,'maxiter':50,'rhobeg':0.1})
x = opt.x
display("Second pass and we have this ... ")
display_result(x)
P.n_calls = 7*50
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints, callback=P.callback,
options={'disp':True,'maxiter':50,'rhobeg':0.01})
display("Third pass and we have this ... ")
display_result(x)
opts[i] = opt
Problem cases¶
The low-end of reject temperature is not coming out well, so I manually determined some initial points.
In [15]:
for T_heat_reject in [303.15]:
print("T_heat_reject = {} K, find max Q".format(T_heat_reject))
tr = T_heat_reject
xB = [400, 1, tr, 3, tr, 5, 285, 4, tr, 0.15]
bdry = system_aqua1.makeBoundary(xB)
P = Problem_2_5_A(bdry)
# For 293 and 298
# x = numpy.array([0.6, 278.0, 306.5, 315.7, 313.0, 371.7])
x = numpy.array([0.50, 278.0, 311.0, 311.2, 307.2, 373.3])
P.n_calls = 7*30
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints, callback=P.callback,
options={'disp':True,'maxiter':50,'rhobeg':0.1})
x = opt.x
display("First pass and we have this ... ")
display_result(x)
P.n_calls = 7*40
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints, callback=P.callback,
options={'disp':True,'maxiter':50,'rhobeg':0.1})
x = opt.x
display("Second pass and we have this ... ")
display_result(x)
P.n_calls = 7*50
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints, callback=P.callback,
options={'disp':True,'maxiter':50,'rhobeg':0.01})
display("Third pass and we have this ... ")
display_result(x)
opts2.at[(tr, 1),'opt'] = dict(opt)
In [14]:
for tr in [303.15]:
for q_frac in [1., 0.95, 0.9, 0.8, 0.5, 0.2, 0.1]:
opts2.at[(tr,q_frac)] = None
Postprocessing - Learning pandas and visualization¶
In [85]:
trendnames = "Q_evap Q_gen m_rich m_refrig UA_gen UA_evap UA_cond UA_rect UA_abs UA_total".split()
#trendtype = numpy.dtype([(name,'f') for name in trendnames])
#trends = numpy.empty_like(t_heat_reject_range,dtype=trendtype)
#trends.fill((None,)*len(trendnames))
trends = pandas.DataFrame(columns=trendnames, index=numpy.arange(len(t_heat_reject_range)))
trends
Out[85]:
In [26]:
for i,(opt, T_heat_reject) in enumerate(zip(opts, t_heat_reject_range)):
print("Case {}, T_heat_reject = {} K".format(i, T_heat_reject))
try:
ch = system_aqua1.makeChiller(opt.x)
trends.Q_evap[i] = ch.Q_evap
trends.Q_gen[i] = ch.Q_gen
trends.m_rich[i] = ch.m_rich
trends.m_refrig[i] = ch.m_refrig
sys = system_aqua1.System(bdry, ch)
trends.UA_gen[i] = sys.df.UA['gen']
trends.UA_evap[i] = sys.df.UA['evap']
except:
pass
In [27]:
trends
Out[27]:
In [28]:
for name in trends:
plt.figure()
trends[name].plot(style='.-')
plt.show()
Step Two: Parameter sweep on output, optimize cost¶
First, let's create a data structure to store the results.
In [22]:
q_fraction_range = [1, 0.999, 0.99, 0.95, 0.9, 0.8, 0.5, 0.2, 0.1]
In [29]:
opts2 = pandas.DataFrame(index=pandas.MultiIndex.from_product([t_heat_reject_range,
q_fraction_range],
names=['T_r','q_frac']),
columns=['Q'] + trendnames + ['opt', ] )
Draw in the results from step one. (TODO: re-arrange all this...)
In [ ]:
for tr,opt in zip(t_heat_reject_range, opts):
opts2.opt[(tr,1)] = opt
ch = system_aqua1.makeChiller(opt.x)
opts2.at[(tr,1),'Q'] = ch.Q_evap
# TODO: Would like to add constraints, eventually
In [17]:
ch = system_aqua1.makeChiller(opt['x'])
opts2.at[(tr,1),'Q'] = ch.Q_evap
In [37]:
#opts2 = opts2.append(pandas.DataFrame(
# index=pandas.MultiIndex.from_product([t_heat_reject_range,
# [0.99, 0.999]],
# names=['T_r','q_frac'])))
#opts2 = opts2.reindex(pandas.MultiIndex.from_product([t_heat_reject_range,
# [1.0, 0.999, 0.99, 0.95, 0.9, 0.8, 0.5, 0.2, 0.1]],
# names=['T_r','q_frac']))
#opts2 = opts2.append(pandas.DataFrame(
# columns=['x0']))
#opts2 = opts2.append(pandas.DataFrame(columns="""deltaT_abs deltaT_cond deltaT_evap deltaT_gen deltaT_rect
# epsilon_abs epsilon_cond epsilon_evap epsilon_gen epsilon_rect""".split()))
In [38]:
opts2
Out[38]:
In [25]:
opts2.iloc[14].opt
Out[25]:
In [116]:
plt.figure(figsize=[8,4])
opts2.Q_evap.plot.bar()
plt.show()
In [70]:
plt.figure(figsize=[8,4])
opts2.UA_total.plot.bar()
plt.ylim(0,200)
plt.show()
In [71]:
fname = 'optimize_trials2_5_data.csv'
if True:
opts2.to_csv(fname)
if False:
opts2 = pandas.read_csv(fname, header=0, index_col=[0,1])
if False:
from numpy import array, inf, nan
# If all the opt entries are __repr__ output from OptimizeResult (oops)
from util import optimize_result_repr_to_dict
for i,o in opts2.opt.items():
if type(o) == str:
print("Making a correction at index {} of opt={}".format(i,o))
if o.startswith('{'):
print("Hm, I think it's a dictionary.")
d = eval(o)
opts2.opt[i] = d
else:
# If all the opt entries are strings as read by read_csv
print("By elimination, probably an OptimizeResult.__repr__() output.")
d = optimize_result_repr_to_dict(o)
opts2.opt[i] = d
In [39]:
a='aardvark'
a.startswith('aa')
Out[39]:
In [47]:
dry_run = True
for i,index in enumerate(opts2.index):
tr,qf = index
q_max = opts2.Q[(tr,1)]
q_goal = q_max * qf
print("Case {}, T_heat_reject = {} K, Q_fraction = {}, Q_goal = {}".format(i, tr, qf, q_goal))
if opts2.opt[index] is numpy.nan or opts2.opt[index] is None:
xB = [400, 1, tr, 3, tr, 5, 285, 4, tr, 0.15]
bdry = system_aqua1.makeBoundary(xB)
P = Problem_2_5_B(bdry,q_goal)
x = numpy.array([0.05, 278.45, tr+7, tr+8, tr+5, 395.15])
if not dry_run:
P.n_calls = 7*30
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints,
options={'disp':True,'maxiter':50,'rhobeg':0.1})
x = opt.x
print("First pass and we have this ... ")
display_result(x)
P.n_calls = 7*40
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints,
options={'disp':True,'maxiter':50,'rhobeg':0.1})
x = opt.x
print("Second pass and we have this ... ")
display_result(x)
P.n_calls = 7*50
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints,
options={'disp':True,'maxiter':50,'rhobeg':0.01})
x = opt.x
print("Third pass and we have this ... ")
display_result(x)
opts2.opt[index] = opt
In [96]:
d=dict()
d[(303.15,0.95)] = array([ 0.47920246, 278.17258004, 311.03722458, 311.29853349,
307.31608332, 373.4217932 ])
d[(303.15,0.9)] = array([ 0.46042825, 278.06261986, 311.12471959, 311.39078344,
307.40903514, 373.42125787])
d[(303.15,0.8)] = array([ 0.40928524, 278.05976775, 311.12458983, 311.38278494,
307.41282053, 373.42343252])
d[(303.15,0.5)] = array([ 2.54650350e-01, 2.78112131e+02, 3.11070930e+02,
3.11363072e+02, 3.07423492e+02, 3.73444409e+02])
d[(303.15,0.2)] = array([ 1.02384404e-01, 2.78082311e+02, 3.11054415e+02,
3.11364635e+02, 3.07485941e+02, 3.73429532e+02])
d[(303.15,0.1)] = array([ 1.00000000e-01, 2.77188460e+02, 3.11580120e+02,
3.11680120e+02, 3.08320712e+02, 3.72906164e+02])
for tr in t_heat_reject_range:
opts2.at[(tr,0.99),'x0'] = opts2.at[(tr,0.95),'opt']['x']
opts2.at[(tr,0.999),'x0'] = opts2.at[(tr,0.95),'opt']['x']
In [99]:
opts2.loc[(293.15,0.99)]
Out[99]:
In [28]:
dry_run = True
for i,index in enumerate(opts2.index):
tr,qf = index
q_max = opts2.Q[(tr,1)]
q_goal = q_max * qf
print("Case {}, T_heat_reject = {} K, Q_fraction = {}, Q_goal = {}".format(i, tr, qf, q_goal))
if opts2.opt[index] in [None, nan]:
xB = [400, 1, tr, 3, tr, 5, 285, 4, tr, 0.15]
bdry = system_aqua1.makeBoundary(xB)
P = Problem_2_5_B(bdry,q_goal)
if index in d:
x = d[index].copy()
else:
x = opts2.opt[(tr,1)]['x'].copy()
print("Initial x = ")
display(x)
if not dry_run:
print("Running first pass optimization ... ")
P.n_calls = 7*10
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints,
options={'disp':True,'maxiter':50,'rhobeg':0.1})
x = opt.x
display_result(x)
print("Running second pass optimization ... ")
P.n_calls = 7*30
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints,
options={'disp':True,'maxiter':50,'rhobeg':0.1})
x = opt.x
display_result(x)
print("Running third pass optimization ... ")
P.n_calls = 7*50
opt = scipy.optimize.minimize(P.objective, x, method="COBYLA",
constraints=P.constraints,
options={'disp':True,'maxiter':50,'rhobeg':0.01})
x = opt.x
display_result(x)
opts2.at[index,'opt'] = dict(opt)
In [104]:
not 1 in [1,2,3]
Out[104]:
Post-processing (again?)¶
In [41]:
# Post-processing
dry_run = True
for i,index in enumerate(opts2.index):
tr,qf = index
q_max = opts2.Q[(tr,1)]
q_goal = q_max * qf
print("Case {}, T_heat_reject = {} K, Q_fraction = {}, Q_goal = {}".format(i, tr, qf, q_goal))
xB = [400, 1, tr, 3, tr, 5, 285, 4, tr, 0.15]
bdry = system_aqua1.makeBoundary(xB)
if not opts2.opt[index] in [None, nan]:
x = opts2.opt[index]['x']
try:
if not dry_run:
ch = system_aqua1.makeChiller(x)
opts2.at[index,'Q_evap'] = ch.Q_evap
opts2.at[index,'Q_gen'] = ch.Q_gen
opts2.at[index,'m_rich'] = ch.m_rich
opts2.at[index,'m_refrig'] = ch.m_refrig
sys = system_aqua1.System(bdry, ch)
opts2.at[index,'UA_gen'] = sys.df.UA['gen']
opts2.at[index,'UA_cond'] = sys.df.UA['cond']
opts2.at[index,'UA_abs'] = sys.df.UA['abs']
opts2.at[index,'UA_rect'] = sys.df.UA['rect']
opts2.at[index,'UA_evap'] = sys.df.UA['evap']
opts2.at[index,'UA_total'] = sys.totalUA
opts2.at[index,'deltaT_gen'] = sys.df.deltaT['gen']
opts2.at[index,'deltaT_cond'] = sys.df.deltaT['cond']
opts2.at[index,'deltaT_abs'] = sys.df.deltaT['abs']
opts2.at[index,'deltaT_rect'] = sys.df.deltaT['rect']
opts2.at[index,'deltaT_evap'] = sys.df.deltaT['evap']
opts2.at[index,'epsilon_gen'] = sys.df.epsilon['gen']
opts2.at[index,'epsilon_cond'] = sys.df.epsilon['cond']
opts2.at[index,'epsilon_abs'] = sys.df.epsilon['abs']
opts2.at[index,'epsilon_rect'] = sys.df.epsilon['rect']
opts2.at[index,'epsilon_evap'] = sys.df.epsilon['evap']
except KeyboardInterrupt:
raise
except:
pass
Results and plots and things¶
In [42]:
display(t_heat_reject_range,q_fraction_range, opts2.head())
Q vs UA, for each heat reject temperature¶
In [119]:
plt.figure()
plt.xlabel('ua total (kW/K)')
plt.ylabel('q cooling (kW)')
for t in t_heat_reject_range:
plt.plot(opts2.UA_total[t], opts2.Q_evap[t], '.-', label="tr={} K".format(t))
plt.legend()
plt.show()
Cooling vs Heat Reject Temperature, for each fraction of max¶
This doesn't have much meaning since the series are relative.
In [120]:
# Group by second index
opts3 = opts2.swaplevel()
In [121]:
plt.figure()
plt.xlabel('t heat reject')
plt.ylabel('q cooling (kW)')
for q_frac in q_fraction_range:
plt.plot(t_heat_reject_range, opts3.Q_evap[q_frac], '.-', label="q_frac={}".format(q_frac))
plt.legend()
plt.show()
Spline things up a bit¶
In [155]:
import scipy.interpolate
x,y,z = (array([tr for (tr,q_frac) in opts2.index]),
array(opts2.UA_total),
array(opts2.Q_evap))
X,Y=numpy.meshgrid(numpy.linspace(293,333),numpy.linspace(0,200))
mask = numpy.logical_not(numpy.isnan(y))
x = x[mask]
y = y[mask]
z = z[mask]
#q_of_t_ua = scipy.interpolate.SmoothBivariateSpline(x,y,z)
#q_of_t_ua = scipy.interpolate.interp2d(x,y,z)
#q_of_t_ua = scipy.interpolate.LSQBivariateSpline(x,y,z,[1,1,1,1,1],[1,1,1,1,1])
#Z=q_of_t_ua(X,Y)
#spline = scipy.interpolate.bisplrep(x,y,z)
Z = scipy.interpolate.griddata(numpy.array([x,y]).T,z,(X,Y))
Z2 = Z.copy()
Z2[numpy.isnan(Z)] = 0
plt.figure()
plt.xlabel('Heat rejection temperature (K)')
plt.ylabel('UA (kW/K)')
plt.title('Q cooling (kW)')
plt.contourf(X,Y,Z)
plt.colorbar()
plt.show()
3D plot of (T,UA,Q) with full extent of data¶
In [122]:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_xlabel('T heat reject (K)')
ax.set_ylabel('ua total (kW/K)')
ax.set_zlabel('q cooling (kW)')
for q_frac in q_fraction_range:
ax.plot(t_heat_reject_range,
array(opts3.UA_total[q_frac]),
array(opts3.Q_evap[q_frac]),
'.-')
for t in t_heat_reject_range:
x = t * numpy.ones_like(q_fraction_range)
ax.plot(x,
array(opts2.UA_total[t]),
array(opts2.Q_evap[t]))
#plt.ylim(0,200)
plt.show()
Cropped 3D wireframe plot¶
In [188]:
# Zoom to relevant data, say, UA < 200
# 3d plot does not handle axis limits well, so instead
# we just need to mask off that data.
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_xlabel('T heat reject (K)')
ax.set_ylabel('ua total (kW/K)')
ax.set_zlabel('q cooling (kW)')
#ax.plot(t_heat_reject_range,
# array(opts3.Q_evap[1.0]),
# zs=200,
# zdir='y')
#for q_frac in q_fraction_range[1:]:
for q_frac in q_fraction_range:
x,y,z = (t_heat_reject_range.copy(),
array(opts3.UA_total[q_frac]),
array(opts3.Q_evap[q_frac]))
mask = y > 200.
#x[mask] = nan
y[mask] = 200
#z[mask] = nan
ax.plot(x,y,z,'b.-')
#ax.plot(array(opts2.UA_total[293.15])[1:],
# array(opts2.Q_evap[293.15])[1:],
# zs=293,
# zdir='x')
for t in t_heat_reject_range:
x,y,z = (t * numpy.ones_like(q_fraction_range),
array(opts2.UA_total[t]),
array(opts2.Q_evap[t]))
mask = y > 200.
#x[mask] = nan
y[mask] = 200
#z[mask] = nan
ax.plot(x,y,z,'r')
ax.set_xlim(290,330)
for i in range(0,X.shape[0],5):
ax.plot(X[i,1:],Y[i,1:],Z2[i,1:],'gray',alpha=0.3)
#ax.contour(X,Y,Z)
plt.show()
In [186]:
# If you got the griddata above ...
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_xlabel('T heat reject (K)')
ax.set_ylabel('ua total (kW/K)')
ax.set_zlabel('q cooling (kW)')
#for i in range(X.shape[0]):
# ax.plot(X[i,:],Y[i,:],Z[i,:],'g',alpha=0.3)
#for j in range(X.shape[1]):
# ax.plot(X[:,j],Y[:,j],Z[:,j],'y',alpha=0.3)
ax.plot_wireframe(X[:,1:],Y[:,1:],Z2[:,1:],color='gray',rstride=4,cstride=4,alpha=0.5)
#ax.plot_surface(X,Y,Z2,cmap=matplotlib.cm.Blues)
ax.contour(X,Y,Z,15)
ax.set_zlim(0,140)
plt.show()
Bubble/color plot¶
In [125]:
# Zoom to relevant data, say, set any UA over threshold to threshold (200)
fig = plt.figure()
plt.xlabel('T heat reject (K)')
plt.ylabel('q cooling (kW)')
plt.title('ua total (kW/K)')
x,y,z = ([tr for (tr,q_frac) in opts2.index],
array(opts2.Q_evap),
array(opts2.UA_total))
plt.scatter(x,y,c=z,vmin=0,vmax=100,cmap='rainbow',s=z,alpha=0.8)
plt.plot(t_heat_reject_range, opts3.Q_evap[1.0])
plt.ylim(0,200)
plt.show()
import scipy.interpolate
In [126]:
# Zoom to relevant data, say, set any UA over threshold to threshold (200)
matplotlib.cm
fig = plt.figure()
plt.xlabel('T heat reject (K)')
plt.ylabel('ua total (kW/K)')
plt.title('q cooling (kW)')
x,y,z = ([tr for (tr,q_frac) in opts2.index],
array(opts2.UA_total),
array(opts2.Q_evap))
plt.scatter(x,y,c=z,cmap='rainbow',s=z,alpha=0.8)
plt.ylim(0,200)
plt.show()
Data viewer¶
In [80]:
%matplotlib inline
import widgetsnbextension
from ipywidgets import interactive
modes = opts2.columns
def f(mode='Q_evap'):
# ch = system_aqua1.makeChiller(xC)
display(mode)
import scipy.interpolate
x,y,z = (array([tr for (tr,q_frac) in opts2.index]),
array(opts2.UA_total),
array(opts2[mode]))
X,Y=numpy.meshgrid(numpy.linspace(293,333),numpy.linspace(0,200))
mask = numpy.logical_not(numpy.isnan(y))
x = x[mask]
y = y[mask]
z = z[mask]
#q_of_t_ua = scipy.interpolate.SmoothBivariateSpline(x,y,z)
#q_of_t_ua = scipy.interpolate.interp2d(x,y,z)
#q_of_t_ua = scipy.interpolate.LSQBivariateSpline(x,y,z,[1,1,1,1,1],[1,1,1,1,1])
#Z=q_of_t_ua(X,Y)
#spline = scipy.interpolate.bisplrep(x,y,z)
Z = scipy.interpolate.griddata(numpy.array([x,y]).T,z,(X,Y))
Z2 = Z.copy()
Z2[numpy.isnan(Z)] = 0
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_xlabel('T heat reject (K)')
ax.set_ylabel('ua total (kW/K)')
ax.set_zlabel(mode)
#for i in range(X.shape[0]):
# ax.plot(X[i,:],Y[i,:],Z[i,:],'g',alpha=0.3)
#for j in range(X.shape[1]):
# ax.plot(X[:,j],Y[:,j],Z[:,j],'y',alpha=0.3)
ax.plot_wireframe(X[:,1:],Y[:,1:],Z2[:,1:],color='gray',rstride=4,cstride=4,alpha=0.5)
#ax.plot_surface(X,Y,Z2,cmap=matplotlib.cm.Blues)
ax.contour(X,Y,Z,15)
#ax.set_zlim(0,140)
plt.show()
interactive(f,
mode=modes)
In [50]:
%%html
<textarea id="TOC-markdown" style="font-family:monospace;width:80%;height:20em;">TOC will be here</textarea>
<script>
$("#TOC-markdown").html(
$('h1,h2,h3,h4').filter(":has(a)").map(function(){
return " ".repeat($(this).prop("tagName")[1])
+ "- <a href='" + encodeURI($(this).children().attr("href")) + "'>"
+ $(this).text() + "</a>";}).get().join("\n")
);
</script>
In [ ]: