1#  ___________________________________________________________________________
2#
3#  Pyomo: Python Optimization Modeling Objects
4#  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC
5#  Under the terms of Contract DE-NA0003525 with National Technology and
6#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
7#  rights in this software.
8#  This software is distributed under the 3-clause BSD License.
9#  ___________________________________________________________________________
10#
11# Model Reference: "A Dynamic HIV-Transmission Model for Evaluating the Costs
12#	and Benefits of Vaccine Programs", D.M. Edwards, R.D. Shachter,
13#	and D.K. Owen 1998, Interfaces
14#
15
16from __future__ import division
17from pyomo.environ import (ConcreteModel, Param, Var, Objective,
18                           Constraint, Set, Expression, Suffix,
19                           value, exp, TransformationFactory)
20from pyomo.dae import ContinuousSet, DerivativeVar
21from pyomo.dae.simulator import Simulator
22from pyomo.contrib.sensitivity_toolbox.sens import sensitivity_calculation
23
24
25def create_model():
26
27    m = ConcreteModel()
28
29    m.tf = Param(initialize=20)
30    m.t = ContinuousSet(bounds=(0,m.tf))
31    m.i = Set(initialize=[0,1,2,3,4,5],ordered=True)
32    m.j = Set(initialize=[0,1],ordered=True)
33    m.ij = Set(initialize=[(0,0),(0,1),(1,0),(1,1),(2,0),(2,1),(3,0),(4,0)],
34                          ordered=True)
35
36    #Set Parameters
37    m.eps = Param(initialize = 0.75, mutable=True)
38
39    m.sig = Param(initialize = 0.15, mutable=True)
40    m.xi  = Param(initialize = 0.983, mutable=True)
41    m.omeg  = Param(initialize = 1.0/10, mutable=True)
42
43    m.Y0   = Param(initialize = 55816.0, mutable=True)
44    m.phi0 = Param(initialize = 0.493, mutable=True)
45
46    d={}
47    d[(0,0)] = 0.0222
48    d[(0,1)] = 0.0222
49    d[(1,0)] = 1/7.1
50    d[(1,1)] = 1/7.1
51    d[(2,0)] = 1/8.1
52    d[(2,1)] = 1/(8.1+5)
53    d[(3,0)] = 1/2.7
54    d[(4,0)] = 1/2.1
55    m.dd = Param(m.ij, initialize=d, mutable=True)
56
57    d_inv={}
58    d_inv[(1,0)] = 7.1
59    d_inv[(2,0)] = 8.1
60    d_inv[(3,0)] = 2.7
61    d_inv[(4,0)] = 2.1
62    m.dd_inv = Param(m.ij, initialize=d_inv, default=0, mutable=True)
63
64    I={}
65    I[(0,0)] =  0.9*m.dd[(0,0)]*m.Y0
66    I[(1,0)] = 0.04*m.dd[(0,0)]*m.Y0
67    I[(2,0)] = 0.04*m.dd[(0,0)]*m.Y0
68    I[(3,0)] = 0.02*m.dd[(0,0)]*m.Y0
69    m.II=Param(m.ij, initialize=I, default=0, mutable=True)
70
71    p={}
72    p[(4,0)] = 0.667
73    m.pp = Param(m.ij, initialize=p, default=2, mutable=True)
74
75    b={}
76    b[(1,0)] = 0.066
77    b[(1,1)] = 0.066
78    b[(2,0)] = 0.066
79    b[(2,1)] = 0.066*(1-0.25)
80    b[(3,0)] = 0.147
81    b[(4,0)] = 0.147
82    m.bb = Param(m.ij,initialize=b, default=0, mutable=True)
83
84    eta00={}
85    eta00[(0,0)] = 0.505
86    eta00[(0,1)] = 0.505
87    eta00[(1,0)] = 0.505
88    eta00[(1,1)] = 0.505
89    eta00[(2,0)] = 0.307
90    eta00[(2,1)] = 0.4803
91    eta00[(3,0)] = 0.235
92    eta00[(4,0)] = 0.235
93    m.eta00 = Param(m.ij, initialize=eta00, mutable=True)
94
95    eta01={}
96    eta01[(0,0)] = 0.505
97    eta01[(0,1)] = 0.6287
98    eta01[(1,0)] = 0.505
99    eta01[(1,1)] = 0.6287
100    eta01[(2,0)] = 0.307
101    eta01[(2,1)] = 0.4803
102    eta01[(3,0)] = 0.235
103    eta01[(4,0)] = 0.235
104    m.eta01 = Param(m.ij, initialize=eta01, mutable=True)
105
106    m.kp  = Param(initialize = 1000.0, mutable=True)
107    m.kt  = Param(initialize = 1000.0, mutable=True)
108    m.rr  = Param(initialize = 0.05, mutable=True)
109
110    c={}
111    c[(0,0)] = 3307
112    c[(0,1)] = 3307
113    c[(1,0)] = 5467
114    c[(1,1)] = 5467
115    c[(2,0)] = 5467
116    c[(2,1)] = 5467
117    c[(3,0)] = 12586
118    c[(4,0)] = 35394
119    m.cc = Param(m.ij, initialize=c, mutable=True)
120
121    q={}
122    q[(0,0)] = 1
123    q[(0,1)] = 1
124    q[(1,0)] = 1
125    q[(1,1)] = 1
126    q[(2,0)] = 0.83
127    q[(2,1)] = 0.83
128    q[(3,0)] = 0.42
129    q[(4,0)] = 0.17
130    m.qq = Param(m.ij, initialize=q, mutable=True)
131
132    m.aa = Param(initialize = 0.0001, mutable=True)
133
134    #Set Variables
135    m.yy = Var(m.t,m.ij)
136    m.L = Var(m.t)
137
138    m.vp = Var(m.t, initialize=0.75, bounds=(0,0.75))
139    m.vt = Var(m.t, initialize=0.75, bounds=(0,0.75))
140
141    m.dyy = DerivativeVar(m.yy, wrt=m.t)
142    m.dL = DerivativeVar(m.L, wrt=m.t)
143
144    def CostFunc(m):
145        return m.L[m.tf]
146    m.cf = Objective(rule=CostFunc)
147
148
149    def _initDistConst(m):
150        return (m.phi0*m.Y0)/sum(m.dd_inv[kk] for kk in m.ij)
151    m.idc = Expression(rule=_initDistConst)
152
153    m.yy[0,(0,0)].fix(value((1-m.phi0)*m.Y0))
154    m.yy[0,(0,1)].fix(0)
155    m.yy[0,(1,0)].fix(value(m.dd_inv[(1,0)]*m.idc))
156    m.yy[0,(1,1)].fix(0)
157    m.yy[0,(2,0)].fix(value(m.dd_inv[(2,0)]*m.idc))
158    m.yy[0,(2,1)].fix(0)
159    m.yy[0,(3,0)].fix(value(m.dd_inv[(3,0)]*m.idc))
160    m.yy[0,(4,0)].fix(value(m.dd_inv[(4,0)]*m.idc))
161    m.L[0].fix(0)
162
163
164    #ODEs
165    def _yy00(m, t):
166        return sum(m.pp[kk]*m.yy[t,kk] for kk in m.ij)*m.dyy[t,(0,0)] == \
167    	   	sum(m.pp[kk]*m.yy[t,kk] for kk in m.ij)*(m.II[(0,0)]-\
168    		(m.vp[t]+m.dd[(0,0)])*m.yy[t,(0,0)]+m.omeg*m.yy[t,(0,1)])-\
169               	m.pp[(0,0)]*sum(m.bb[kk]*m.eta00[kk]*m.pp[kk]*m.yy[t,kk]
170                                for kk in m.ij)*m.yy[t,(0,0)]
171    m.yy00DiffCon = Constraint(m.t, rule=_yy00)
172
173    def _yy01(m, t):
174        return sum(m.pp[kk]*m.yy[t,kk] for kk in m.ij)*m.dyy[t,(0,1)] == \
175               	sum(m.pp[kk]*m.yy[t,kk]
176                    for kk in m.ij)*(m.vp[t]*m.yy[t,(0,0)]-
177                                    (m.dd[(0,0)]+m.omeg)*m.yy[t,(0,1)])-\
178                m.pp[(0,1)]*(1-m.eps)*sum(m.bb[kk]*m.eta01[kk]*
179                                          m.pp[kk]*m.yy[t,kk]
180                                          for kk in m.ij)*m.yy[t,(0,1)]
181    m.yy01DiffCon = Constraint(m.t, rule=_yy01)
182
183    def _yy10(m, t):
184        return sum(m.pp[kk]*m.yy[t,kk] for kk in m.ij)*m.dyy[t,(1,0)] == \
185               	sum(m.pp[kk]*m.yy[t,kk] for kk in m.ij)*\
186                (m.II[(1,0)]-((m.sig*m.xi)+m.vp[t]+m.dd[(1,0)]+m.dd[(0,0)])*
187                  m.yy[t,(1,0)]+m.omeg*m.yy[t,(1,1)]
188                )+m.pp[(0,0)]*sum(m.bb[kk]*m.eta00[kk]*
189                                  m.pp[kk]*m.yy[t,kk]
190                                  for kk in m.ij)*m.yy[t,(0,0)]
191    m.yy10DiffCon = Constraint(m.t, rule=_yy10)
192
193    def _yy11(m, t):
194        return sum(m.pp[kk]*m.yy[t,kk] for kk in m.ij)*m.dyy[t,(1,1)] == \
195               	sum(m.pp[kk]*m.yy[t,kk] for kk in m.ij)*(m.vp[t]*m.yy[t,(1,0)]-\
196                (m.omeg+(m.sig*m.xi)+m.dd[(1,1)]+m.dd[(0,0)])*m.yy[t,(1,1)])+\
197    		m.pp[(0,1)]*(1-m.eps)*sum(m.bb[kk]*m.eta01[kk]*
198                                          m.pp[kk]*m.yy[t,kk]
199                                          for kk in m.ij)*m.yy[t,(0,1)]
200    m.yy11DiffCon = Constraint(m.t, rule=_yy11)
201
202    def _yy20(m, t):
203        return m.dyy[t,(2,0)] == \
204    		m.II[(2,0)]+m.sig*m.xi*(m.yy[t,(1,0)]+m.yy[t,(1,1)])-\
205    		(m.vt[t]+m.dd[(2,0)]+m.dd[(0,0)])*m.yy[t,(2,0)]
206    m.yy20DiffCon = Constraint(m.t, rule=_yy20)
207
208    def _yy21(m, t):
209        return m.dyy[t,(2,1)] == \
210    		m.vt[t]*m.yy[t,(2,0)]-(m.dd[(2,1)]+m.dd[(0,0)])*m.yy[t,(2,1)]
211    m.yy21DiffCon = Constraint(m.t, rule=_yy21)
212
213    def _yy30(m, t):
214        return m.dyy[t,(3,0)] == \
215    		m.II[(3,0)]+m.dd[(1,0)]*m.yy[t,(1,0)]+\
216                m.dd[(1,1)]*m.yy[t,(1,1)]+m.dd[(2,0)]*m.yy[t,(2,0)]+\
217                m.dd[(2,1)]*m.yy[t,(2,1)]-\
218                (m.dd[(3,0)]+m.dd[(0,0)])*m.yy[t,(3,0)]
219    m.yy30DiffCon = Constraint(m.t, rule=_yy30)
220
221    def _yy40(m, t):
222        return m.dyy[t, (4,0)] == \
223    		m.dd[(3,0)]*m.yy[t,(3,0)]-(m.dd[(4,0)]+\
224                m.dd[(0,0)])*m.yy[t,(4,0)]
225    m.yy40DiffCon = Constraint(m.t, rule=_yy40)
226
227    def _L(m, t):
228        return m.dL[t] == \
229            exp(-m.rr*t)*(m.aa*(m.kp*m.vp[t]*(m.yy[t,(0,0)]+m.yy[t,(1,0)]) \
230            +(m.kt*m.vt[t]*m.yy[t,(2,0)])+sum(m.cc[kk]*m.yy[t,kk]
231                                              for kk in m.ij)) \
232            -(1-m.aa)*sum(m.qq[kk]*m.yy[t,kk] for kk in m.ij))
233    m.LDiffCon = Constraint(m.t, rule=_L)
234
235    return m
236
237
238def initialize_model(m,n_sim,n_nfe,n_ncp):
239    vp_profile = {0:0.75}
240    vt_profile = {0:0.75}
241
242
243    m.u_input = Suffix(direction=Suffix.LOCAL)
244    m.u_input[m.vp] = vp_profile
245    m.u_input[m.vt] = vt_profile
246
247    sim = Simulator(m, package='scipy')
248    tsim, profiles = sim.simulate(numpoints=n_sim, varying_inputs=m.u_input)
249
250
251    discretizer = TransformationFactory('dae.collocation')
252    discretizer.apply_to(m,nfe=n_nfe,ncp=n_ncp,scheme='LAGRANGE-RADAU')
253
254    sim.initialize_model()
255
256
257
258if __name__=='__main__':
259    m = create_model()
260    initialize_model(m,10,5,1)
261
262    m.epsDelta = Param(initialize = 0.75001)
263
264    q_del={}
265    q_del[(0,0)] = 1.001
266    q_del[(0,1)] = 1.002
267    q_del[(1,0)] = 1.003
268    q_del[(1,1)] = 1.004
269    q_del[(2,0)] = 0.83001
270    q_del[(2,1)] = 0.83002
271    q_del[(3,0)] = 0.42001
272    q_del[(4,0)] = 0.17001
273    m.qqDelta = Param(m.ij, initialize=q_del)
274
275    m.aaDelta = Param(initialize = .0001001)
276
277    m_sipopt = sensitivity_calculation('sipopt', m,[m.eps,m.qq,m.aa],
278                        [m.epsDelta,m.qqDelta,m.aaDelta],
279                        tee = True)
280