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