-
Notifications
You must be signed in to change notification settings - Fork 0
/
plat_Cfswitch.py
168 lines (143 loc) · 4.55 KB
/
plat_Cfswitch.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
"""
This script calculates the Pareto front for the Cf-switch mechanism using
multi-objective optimzation methods provided in the package Platypus
(https://platypus.readthedocs.io/en/latest/).
Written by Wylie Stroberg in 2018
"""
import matplotlib as mpl
mpl.use('TKAgg')
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import simps
import time
from platypus import NSGAII, NSGAIII, SPEA2, Problem, Real, ProcessPoolEvaluator
from delayed_nf_upr_piecewise_cf_switch import *
#------------------------------------------------------------
#------------------------------------------------------------
#------------------------------------------------------------
def objective_function(dde,params):
sol = run_dde(dde,params)
c = sol['c']
u = sol['u']
t = sol['t']
USS, CSS = calc_steady_state(params)
Cin = c[-1]-c[0] + simps(c,t)
Udiff = u-USS
Uexcess = np.where(Udiff<0,0,Udiff)
Utot = simps(Uexcess,t)
# Normalize
cin = Cin/(t[-1])
utot = Utot/(USS*t[-1]) + 1
return [cin,utot]
#------------------------------------------------------------
#------------------------------------------------------------
if __name__=="__main__":
# Physical Parameters
KCAT = 8.15e-4
KM = 1.1e4
VU = 200.
VC = 60.
B = 1.85e-4
G = 5.0
TAU_UPR = 15.*60.
kcat = KCAT+B
# Non-dimensional parameters
nu = VU/VC
alpha = kcat/B
beta = KM*B/VC
tau_upr = TAU_UPR*B
# Pulse Parameters
pulse_type = 'med'
if pulse_type=='chronic':
d = 1.0
tau_p = 1000.0*60.*B
elif pulse_type=='acute':
d = 100.0
tau_p = 10.0*60.*B
elif pulse_type=='med':
d = 1.0
tau_p = 50.0*60.*B
totaluin = 0.5
d = 1.0e0
tau_p = totaluin/d
teq = 30.*tau_upr
# Response function variables (values not used)
m = 1.0
cmax = 0.6
# Initialize dde solver
params = {
'alpha' : alpha,
'beta' : beta,
'nu' : nu,
'G' : G,
'cmax' : cmax,
'm' : m,
'tau_upr' : tau_upr,
'd' : d,
'tau_p' : tau_p,
'teq' : teq
}
#tfinal = 2000.*60 + teq
tfinal = teq + tau_p + 30.*tau_upr
dtmax = 0.003
dde = initialize_dde_nf_upr(params=params,tfinal=tfinal,dtmax=dtmax)
USS, CSS = calc_steady_state(params)
CF_SS = CSS*(1.-USS/(beta+USS))
# Define objective function w/ constraints
def uswitch_objective(vars):
params['cmax'] = vars[0]
params['m'] = vars[1]
obj = objective_function(dde,params)
pad = 0.001
constraint = vars[0] - CF_SS*(1.+pad)
return (obj,[constraint])
#test_vars = (26.08718718667782, 577.8399829310629)
#print(uswitch_objective(test_vars))
# Define multi-objective optimization problem object
cmax_min = CF_SS*1e-3
cmax_max = CF_SS*0.99 # for short pulse
mmin = 0.01
mmax = 1000.0
problem = Problem(2,2,1)
problem.types[:] = [Real(cmax_min,cmax_max),Real(mmin,mmax)]
problem.constraints[:] = "<=0"
problem.function = uswitch_objective
# Run optimization in parallel
n_proc = 4
n_evals = 100000
alg_name = 'NSGAII'
#alg_name = 'NSGAIII'
#alg_name = 'SPEA2'
start_time = time.time()
with ProcessPoolEvaluator(n_proc) as evaluator:
if alg_name=='NSGAII':
algorithm = NSGAII(problem,evaluator=evaluator)
elif alg_name=='NSGAIII':
divs = 50
alg_name=alg_name+str(divs)
algorithm = NSGAIII(problem,divisions_outer=divs,evaluator=evaluator)
elif alg_name=='SPEA2':
algorithm = SPEA2(problem,evaluator=evaluator)
else:
print("Using NSGAII algorithm by default")
algorithm = NSGAII(problem,evaluator=evaluator)
algorithm.run(n_evals)
print("Total run time = {:.2f} hours".format((time.time()-start_time)/3600.))
# Save Results
save_dir = "./Data/Cf_switch/plat/"
pareto_file_name = save_dir + "pareto_front_{0}_{1}_{2}_".format(d,tau_p,n_evals) + alg_name + "_dtmax{1}.npy".format(n_proc,dtmax)
points = []
costs = []
for s in algorithm.result:
points.append(s.variables[:])
costs.append(s.objectives[:])
points = np.array(points)
costs = np.array(costs)
pareto_front = np.hstack((points,costs))
np.save(pareto_file_name,pareto_front)
# Plot Pareto set
plt.scatter([s.objectives[0] for s in algorithm.result],
[s.objectives[1] for s in algorithm.result])
plt.xlabel("$f_1(x)$")
plt.ylabel("$f_2(x)$")
plt.show()