From ff7f0ee656091b28c3cfb4adf577dd9c6224a630 Mon Sep 17 00:00:00 2001
From: Gabriel Patron <g2patron@uwaterloo.ca>
Date: Mon, 12 Aug 2024 05:19:52 -0400
Subject: [PATCH] Upload New File
---
MPC_RAS_tempdist_git.py | 1995 +++++++++++++++++++++++++++++++++++++++
1 file changed, 1995 insertions(+)
create mode 100644 MPC_RAS_tempdist_git.py
diff --git a/MPC_RAS_tempdist_git.py b/MPC_RAS_tempdist_git.py
new file mode 100644
index 0000000..dd70fe6
--- /dev/null
+++ b/MPC_RAS_tempdist_git.py
@@ -0,0 +1,1995 @@
+
+# everything has been checked
+
+from pyomo.environ import *
+from pyomo.dae import *
+from pyomo.opt import SolverFactory
+import numpy as mp
+import matplotlib.pyplot as plt
+import math
+InitGUStates={0: 0.0004941262758492059, 1: 0.0026610418262743765, 2: 0.0005953489739768516, 3: 0.007999999921084513, 4: 0.008, 5: 0, 6: 0.0012704702721796658, 7: 0.0006551514887411368, 8: 0.024912867884266228, 9: 0.0004521422131039016, 10: 0.01048087596374192, 11: 0.0033743585516853034, 12: 0.021094991040364303, 13: 3.258529553051713e-05, 14: 0.2548137888637069, 15: 0.08203828436634969, 16: 0.002858920928472409, 17: 0.00032623896911251616, 18: 0.006184272269933473, 19: 0.019571519207748106, 20: 0.0004941262758492059, 21: 0.0026610418262743765, 22: 0.0005953489739768516, 23: 0.007999999921084513, 24: 0.008, 25: 0, 26: 0.0012704702721796658, 27: 0.0006551514887411368, 28: 0.024912867884266228, 29: 0.0004521422131039016, 30: 0.01048087596374192, 31: 0.0033743585516853034, 32: 0.021094991040364303, 33: 3.258529553051713e-05, 34: 0.2548137888637069, 35: 0.08203828436634969, 36: 0.002858920928472409, 37: 0.00032623896911251616, 38: 0.006184272269933473, 39: 0.019571519207748106, 40: 0.0006248485572695803, 41: 0.002697183223242683, 42: 0.0006262134488224627, 43: 0.007999999922041812, 44: 0.008, 45: 0.001476084484174653, 46: 0.0012976268154275138, 47: 0.024912867884244745, 48: 0.0008619423141437846, 49: 0, 50: 3.2301907373349344e-05, 51: 0.0033450122867474297, 52: 0.010389725435038574, 53: 0.020911531222781592, 54: 0.04707215999973613, 55: 0.019415519999952287, 56: 0.009498240000433996, 57: 0.06460622112977724, 58: 0.19318526906403832, 59: 0.06460622112967648, 60: 0.12351189333605922, 61: 0.0002822708778531383, 62: 0.0002822708778531383, 63: 0.0002798160125124163, 64: 0.00047569827657743363, 65: 0.00047569827657743363, 66: 0.0004715612068880054, 67: 0.011565295475379671, 68: 0.011565295475379671, 69: 45, 70: 2300}
+
+
+def MPC_FISH(h,P,States,inputs,n_coll,flag,xsp,ghaza):
+ m=ConcreteModel()
+ m.t=ContinuousSet(bounds=(0,h*P))
+
+ m.F=Var(m.t,initialize=ghaza[1],bounds=(0,3)) #manipulated variable
+ m.kl=Var(m.t,bounds=(0,1),initialize=inputs[1]) #manipulated variable
+ m.klBR1=Var(m.t,bounds=(0,1),initialize=inputs[2]) #manipulated variable
+ m.klBR2=Var(m.t,bounds=(0,1),initialize=inputs[4]) #manipulated variable
+ m.Q4=Var(m.t,bounds=(0,3),initialize=inputs[3]) #manipulated variable
+
+ m.n0=Param(initialize=States[70],mutable = True) # Number of fish when the fish size reaches 6
+ m.agro=Param(initialize=0.4)
+ m.hgro=Param(initialize= 0.08)
+ m.T=Param(initialize=19.5, mutable = True) # Temperature (C)
+ m.Tmingro=Param(initialize=1)
+ m.Tmaxgro=Param(initialize=25)
+ m.Toptgro=Param(initialize=16)
+ m.ratioo=Param(initialize=1)
+ m.Q =Param(initialize=172.8 ) # 2*172.8 or the flow from the MBBs to the fish tank
+ m.Q1 =Param(initialize=345.6)
+
+ def _tau(m):
+ return exp(-4.6*((m.T-m.Toptgro)/(m.Tmaxgro-m.Toptgro))**4)
+ m.tau=Expression(rule=_tau)
+
+ m.jgro=Param(initialize=0.0132)
+ m.Kmingro=Param(initialize=0.00133)
+ m.mgro=Param(initialize=0.67)
+ m.ngro=Param(initialize=0.81)
+ m.Kgro=m.Kmingro*exp(m.jgro*(m.T-m.Tmingro))
+
+
+ m.DOcritgro=Param(initialize=6.98) #DOmaxFI
+ m.DOmingro=Param(initialize=5.3 ) #LOS
+
+ m.Acritgro=Param(initialize=0.0125)
+ m.Amaxgro=Param(initialize=1.4)
+
+ m.fgro=Param(initialize=1)
+
+ m.tav1=Param(initialize=0.211) #Feed residence time(day)
+ m.V=Param(initialize=5.5) # Fish tank volume (m3)
+ m.IBW=Param(initialize=45) # initial body weight of fish(g)
+
+ ### Feed composition
+ m.prot=Param(initialize=0.45)
+ m.lip=Param(initialize=0.24)
+ m.NFE=Param(initialize=0.154)
+ m.fiber=Param(initialize=0.027)
+
+ m.ADCprot=Param(initialize=0.94)
+ m.ADClip=Param(initialize=0.91)
+ m.ADCNFE=Param(initialize=0.66)
+
+
+ m.muh =Param(initialize=4)
+ m.bH=Param(initialize=0.3)
+ m.muAOB =Param(initialize=0.63)
+ m.muNOB=Param(initialize=1.04)
+ m.bAOB =Param(initialize=0.05)
+ m.bNOB =Param(initialize=0.05)
+ m.kh =Param(initialize=3)
+ m.Kx=Param(initialize=0.1)
+
+ m.etano2=Param(initialize=0.8)
+ m.etano3=Param(initialize=0.8)
+
+ # Kinetic paramters
+ m.KS=Param(initialize=10*(1e-3))
+ m.KOH =Param(initialize=0.2*(1e-3),mutable=True)
+ m.KNOx=Param(initialize=0.5* (1e-3))
+ m.KNO2=Param(initialize=1* (1e-3))
+ m.etag =Param(initialize=0.8)
+
+ m.Ka =Param(initialize=0.05*1000)
+ m.KNH =Param(initialize=1*(1e-3))
+ m.KNHI =Param(initialize=5*(1e-3))
+
+ m.KOA=Param(initialize=0.4*(1e-3),mutable=True )
+
+
+ # Stoichiometric parameters
+ m.YAOB =Param(initialize=0.21)
+ m.YNOB =Param(initialize=0.03)
+
+ m.YH =Param(initialize=0.67)
+ m.fp=Param(initialize=0.08)
+ m.iXB=Param(initialize=0.08)
+ m.iXP=Param(initialize=0.06)
+
+ # Inlet concentration
+ m.Hint=Param(initialize=0)
+ m.Aint=Param(initialize=0)
+
+
+
+
+
+ # 0D-BFM parameters
+ m.icv =Param(initialize=1)
+ m.a=Param(initialize=0)
+ m.Vw1 =Param(initialize=0.4)
+
+
+ # Decision variables
+ m.kka=Param(initialize=300)
+ m.kkd=Param(initialize=0.9)
+ m.k=Param(initialize=7.48e-5 )
+ m.loss=Param(initialize=0.03)
+
+ def _nu(m,i):
+ return 1
+ m.nu=Expression(m.t,rule=_nu)
+
+
+ m.pH=Param(initialize=7)
+ m.pKa=Param(initialize=9.245)
+
+ m.cgro=Param(initialize=0.7)
+ m.dgro=Param(initialize=0.9)
+ m.epsgro=Param(initialize=0.05)
+ m.Tcgro=Param(initialize=130)
+
+ m.ni0= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni1= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni2= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni3= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni4= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni5= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni6= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni7= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni8= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni9= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni10= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni11= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni12= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni13= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni14= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni15= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni16= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni17= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni18= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni19= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni20= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni21= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni22= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni23= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni24= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni25= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni26= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni27= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni28= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni29= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni30= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni31= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni32= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni33= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni34= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni35= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni36= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni37= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni38= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni39= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni40= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni41= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni42= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni43= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni44= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni45= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni46= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni47= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni48= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni49= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni50= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni51= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni52= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni53= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni54= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni55= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni56= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni57= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni58= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni59= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni60= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni61= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni62= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni63= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni64= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni65= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni66= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni67= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni68= Param(initialize=-1, mutable=True, within=Integers)
+ m.ni69= Param(initialize=-1, mutable=True, within=Integers)
+
+
+ m.z4=Param(initialize=0.008)
+ m.z5=Param(initialize=0.008)
+ m.z24=Param(initialize=0.008)
+ m.z25=Param(initialize=0.008)
+ m.z44=Param(initialize=0.008)
+ m.z49=Param(initialize=0.008)
+
+ m.z0=Var(m.t,initialize=InitGUStates[0],within=NonNegativeReals)
+ m.z1=Var(m.t,initialize=InitGUStates[1],within=NonNegativeReals)
+ m.z2=Var(m.t,initialize=InitGUStates[2],within=NonNegativeReals)
+ m.z3=Var(m.t,initialize=InitGUStates[3],within=NonNegativeReals)
+
+ m.z6=Var(m.t,initialize=InitGUStates[6],within=NonNegativeReals)
+ m.z7=Var(m.t,initialize=InitGUStates[7],within=NonNegativeReals)
+ m.z8=Var(m.t,initialize=InitGUStates[8],within=NonNegativeReals)
+ m.z9=Var(m.t,initialize=InitGUStates[9],within=NonNegativeReals)
+ m.z10=Var(m.t,initialize=InitGUStates[10],within=NonNegativeReals)
+ m.z11=Var(m.t,initialize=InitGUStates[11],within=NonNegativeReals)
+ m.z12=Var(m.t,initialize=InitGUStates[12],within=NonNegativeReals)
+ m.z13=Var(m.t,initialize=InitGUStates[13],within=NonNegativeReals)
+ m.z14=Var(m.t,initialize=InitGUStates[14],within=NonNegativeReals)
+ m.z15=Var(m.t,initialize=InitGUStates[15],within=NonNegativeReals)
+ m.z16=Var(m.t,initialize=InitGUStates[16],within=NonNegativeReals)
+ m.z17=Var(m.t,initialize=InitGUStates[17],within=NonNegativeReals)
+ m.z18=Var(m.t,initialize=InitGUStates[18],within=NonNegativeReals)
+ m.z19=Var(m.t,initialize=InitGUStates[19],within=NonNegativeReals)
+ m.z20=Var(m.t,initialize=InitGUStates[20],within=NonNegativeReals)
+ m.z21=Var(m.t,initialize=InitGUStates[21],within=NonNegativeReals)
+ m.z22=Var(m.t,initialize=InitGUStates[22],within=NonNegativeReals)
+ m.z23=Var(m.t,initialize=InitGUStates[23],within=NonNegativeReals)
+
+ m.z26=Var(m.t,initialize=InitGUStates[26],within=NonNegativeReals)
+ m.z27=Var(m.t,initialize=InitGUStates[27],within=NonNegativeReals)
+ m.z28=Var(m.t,initialize=InitGUStates[28],within=NonNegativeReals)
+ m.z29=Var(m.t,initialize=InitGUStates[29],within=NonNegativeReals)
+ m.z30=Var(m.t,initialize=InitGUStates[30],within=NonNegativeReals)
+ m.z31=Var(m.t,initialize=InitGUStates[31],within=NonNegativeReals)
+ m.z32=Var(m.t,initialize=InitGUStates[32],within=NonNegativeReals)
+ m.z33=Var(m.t,initialize=InitGUStates[33],within=NonNegativeReals)
+ m.z34=Var(m.t,initialize=InitGUStates[34],within=NonNegativeReals)
+ m.z35=Var(m.t,initialize=InitGUStates[35],within=NonNegativeReals)
+ m.z36=Var(m.t,initialize=InitGUStates[36],within=NonNegativeReals)
+ m.z37=Var(m.t,initialize=InitGUStates[37],within=NonNegativeReals)
+ m.z38=Var(m.t,initialize=InitGUStates[38],within=NonNegativeReals)
+ m.z39=Var(m.t,initialize=InitGUStates[39],within=NonNegativeReals)
+
+ m.z40=Var(m.t,initialize=InitGUStates[40],within=NonNegativeReals)
+ m.z41=Var(m.t,initialize=InitGUStates[41],within=NonNegativeReals)
+ m.z42=Var(m.t,initialize=InitGUStates[42],within=NonNegativeReals)
+ m.z43=Var(m.t,initialize=InitGUStates[43],within=NonNegativeReals)
+ m.z45=Var(m.t,initialize=InitGUStates[45],within=NonNegativeReals)
+ m.z46=Var(m.t,initialize=InitGUStates[46],within=NonNegativeReals)
+ m.z47=Var(m.t,initialize=InitGUStates[47],within=NonNegativeReals)
+ m.z48=Var(m.t,initialize=InitGUStates[48],within=NonNegativeReals)
+ m.z50=Var(m.t,initialize=InitGUStates[50],within=NonNegativeReals)
+ m.z51=Var(m.t,initialize=InitGUStates[51],within=NonNegativeReals)
+ m.z52=Var(m.t,initialize=InitGUStates[52],within=NonNegativeReals)
+ m.z53=Var(m.t,initialize=InitGUStates[53],within=NonNegativeReals)
+ m.z54=Var(m.t,initialize=InitGUStates[54],within=NonNegativeReals)
+ m.z55=Var(m.t,initialize=InitGUStates[55],within=NonNegativeReals)
+ m.z56=Var(m.t,initialize=InitGUStates[56],within=NonNegativeReals)
+ m.z57=Var(m.t,initialize=InitGUStates[57],within=NonNegativeReals)
+ m.z58=Var(m.t,initialize=InitGUStates[58],within=NonNegativeReals)
+ m.z59=Var(m.t,initialize=InitGUStates[59],within=NonNegativeReals)
+ m.z60=Var(m.t,initialize=InitGUStates[60],within=NonNegativeReals)
+ m.z61=Var(m.t,initialize=InitGUStates[61],within=NonNegativeReals)
+ m.z62=Var(m.t,initialize=InitGUStates[62],within=NonNegativeReals)
+ m.z63=Var(m.t,initialize=InitGUStates[63],within=NonNegativeReals)
+ m.z64=Var(m.t,initialize=InitGUStates[64],within=NonNegativeReals)
+ m.z65=Var(m.t,initialize=InitGUStates[65],within=NonNegativeReals)
+ m.z66=Var(m.t,initialize=InitGUStates[66],within=NonNegativeReals)
+ m.z67=Var(m.t,initialize=InitGUStates[67],within=NonNegativeReals)
+ m.z68=Var(m.t,initialize=InitGUStates[68],within=NonNegativeReals)
+ m.z69=Var(m.t,initialize=InitGUStates[69],within=NonNegativeReals)
+ m.z70=Var(m.t,initialize=InitGUStates[70],within=NonNegativeReals)
+
+
+ m.dNHdt= DerivativeVar(m.z0, wrt=m.t)
+ m.dSNDdt= DerivativeVar(m.z1, wrt=m.t)
+ m.dXNDdt= DerivativeVar(m.z2, wrt=m.t)
+ m.dz02= DerivativeVar(m.z3, wrt=m.t)
+ m.dzSdt= DerivativeVar(m.z6, wrt=m.t)
+ m.dzXdt= DerivativeVar(m.z7, wrt=m.t)
+ m.dzSIdt= DerivativeVar(m.z8, wrt=m.t)
+ m.dzXIdt= DerivativeVar(m.z9, wrt=m.t)
+ m.dXHdt= DerivativeVar(m.z10, wrt=m.t)
+ m.dAdt= DerivativeVar(m.z11, wrt=m.t)
+ m.dzNOdt= DerivativeVar(m.z12, wrt=m.t)
+ m.dXPdt= DerivativeVar(m.z13, wrt=m.t)
+ m.dXHbdt= DerivativeVar(m.z14, wrt=m.t)
+ m.dAbdt= DerivativeVar(m.z15, wrt=m.t)
+ m.dzXbdt= DerivativeVar(m.z16, wrt=m.t)
+ m.dXNDbdt= DerivativeVar(m.z17, wrt=m.t)
+ m.dzXIbdt= DerivativeVar(m.z18, wrt=m.t)
+ m.dXPbdt= DerivativeVar(m.z19, wrt=m.t)
+
+
+
+ m.dNHdt2= DerivativeVar(m.z20, wrt=m.t)
+ m.dSNDdt2= DerivativeVar(m.z21, wrt=m.t)
+ m.dXNDdt2= DerivativeVar(m.z22, wrt=m.t)
+ m.dz022= DerivativeVar(m.z23, wrt=m.t)
+ m.dzSdt2= DerivativeVar(m.z26, wrt=m.t)
+ m.dzXdt2= DerivativeVar(m.z27, wrt=m.t)
+ m.dzSIdt2= DerivativeVar(m.z28, wrt=m.t)
+ m.dzXIdt2= DerivativeVar(m.z29, wrt=m.t)
+ m.dXHdt2= DerivativeVar(m.z30, wrt=m.t)
+ m.dAdt2= DerivativeVar(m.z31, wrt=m.t)
+ m.dzNOdt2= DerivativeVar(m.z32, wrt=m.t)
+ m.dXPdt2= DerivativeVar(m.z33, wrt=m.t)
+ m.dXHbdt2= DerivativeVar(m.z34, wrt=m.t)
+ m.dAbdt2= DerivativeVar(m.z35, wrt=m.t)
+ m.dzXbdt2= DerivativeVar(m.z36, wrt=m.t)
+ m.dXNDbdt2= DerivativeVar(m.z37, wrt=m.t)
+ m.dzXIbdt2= DerivativeVar(m.z38, wrt=m.t)
+ m.dXPbdt2= DerivativeVar(m.z39, wrt=m.t)
+
+
+
+ m.dy1dt= DerivativeVar(m.z40, wrt=m.t)
+ m.dy2dt= DerivativeVar(m.z41, wrt=m.t)
+ m.dy3dt= DerivativeVar(m.z42, wrt=m.t)
+ m.dy4dt= DerivativeVar(m.z43, wrt=m.t)
+ m.dy6dt= DerivativeVar(m.z45, wrt=m.t)
+ m.dy7dt= DerivativeVar(m.z46, wrt=m.t)
+ m.dy8dt= DerivativeVar(m.z47, wrt=m.t)
+ m.dy9dt= DerivativeVar(m.z48, wrt=m.t)
+ m.dy11dt= DerivativeVar(m.z50, wrt=m.t)
+ m.dy12dt= DerivativeVar(m.z51, wrt=m.t)
+ m.dy13dt= DerivativeVar(m.z52, wrt=m.t)
+ m.dy14dt= DerivativeVar(m.z53, wrt=m.t)
+
+
+ m.dx1dt= DerivativeVar(m.z54, wrt=m.t)
+ m.dx2dt= DerivativeVar(m.z55, wrt=m.t)
+ m.dx3dt= DerivativeVar(m.z56, wrt=m.t)
+ m.dx6dt= DerivativeVar(m.z57, wrt=m.t)
+ m.dx7dt= DerivativeVar(m.z58, wrt=m.t)
+ m.dx8dt= DerivativeVar(m.z59, wrt=m.t)
+ m.dx9dt= DerivativeVar(m.z60, wrt=m.t)
+
+
+
+ m.ddzNOdt= DerivativeVar(m.z61, wrt=m.t)
+ m.ddzNOdt2= DerivativeVar(m.z62, wrt=m.t)
+ m.ddy14dt= DerivativeVar(m.z63, wrt=m.t)
+ m.ddAdt= DerivativeVar(m.z64, wrt=m.t)
+ m.ddAdt2= DerivativeVar(m.z65, wrt=m.t)
+ m.ddy12dt= DerivativeVar(m.z66, wrt=m.t)
+ m.ddAbdt= DerivativeVar(m.z67, wrt=m.t)
+ m.ddAbdt2= DerivativeVar(m.z68, wrt=m.t)
+ m.dWdt= DerivativeVar(m.z69, wrt=m.t)
+
+ def _delta(m,i):
+ return 1/(1+exp(-5*(m.z43[i]*1000-6.14)))
+ m.delta=Expression(m.t,rule=_delta)
+
+ def _n(m,i):
+ return m.z70[i] == m.n0*exp(-(m.z69[i]**(-0.605))*i*0.0037)
+ m.n=Constraint(m.t,rule=_n)
+
+ def _CODsolid(m):
+ return 1.77* m.prot * (1- m.ADCprot) + 2.88* m.lip * (1- m.ADClip) + 1.16* m.NFE * (1- m.ADCNFE) + 1.16 * m.fiber * (1- m.ADCNFE)
+ m.CODsolid=Expression(rule=_CODsolid)
+
+ def _CODsolidactual(m):
+ return 0.9*m.CODsolid
+ m.CODsolidactual=Expression(rule=_CODsolidactual)
+
+ def _CODfeed(m):
+ return 1.77* m.prot + 2.88* m.lip + 1.16* m.NFE + 1.16 * m.fiber
+ m.CODfeed=Expression(rule=_CODfeed)
+
+ def _CODfeedactual(m):
+ return 0.74 * m.CODfeed
+ m.CODfeedactual=Expression(rule=_CODfeedactual)
+
+ def _Ningested(m):
+ return m.prot/6.25
+ m.Ningested=Expression(rule=_Ningested)
+
+ def _biom(m,i):
+ return (m.z69[i]*m.z70[i])/1000
+ m.biom=Expression(m.t,rule=_biom)
+
+ def _DWN1(m,i):
+ return (1-m.loss)*m.F[i]* m.Ningested *0.337
+ m.DWN1=Expression(m.t,rule=_DWN1)
+
+ m.Loss1=Param(initialize=0)
+
+ def _DWN2(m,i):
+ return (1-m.loss)*m.F[i]* m.Ningested *0.139
+ m.DWN2=Expression(m.t,rule=_DWN2)
+
+ def _Loss2(m,i):
+ return m.loss *m.F[i]* m.Ningested *0.29
+ m.Loss2=Expression(m.t,rule=_Loss2)
+
+ def _DWN3(m,i):
+ return (1-m.loss)* m.F[i]* m.Ningested *0.068
+ m.DWN3=Expression(m.t,rule=_DWN3)
+
+ def _Loss3(m,i):
+ return m.loss *m.F[i]* m.Ningested *0.71
+ m.Loss3=Expression(m.t,rule=_Loss3)
+
+ def _DWN6(m,i):
+ return (1-m.loss)* m.F[i]* 0.204 *m.CODsolidactual
+ m.DWN6=Expression(m.t,rule=_DWN6)
+
+ def _Loss6(m,i):
+ return m.loss *m.F[i] *0.145 *m.CODfeedactual
+ m.Loss6=Expression(m.t,rule=_Loss6)
+
+ def _DWN7(m,i):
+ return (1-m.loss)* m.F[i]* 0.61 *m.CODsolidactual
+ m.DWN7=Expression(m.t,rule=_DWN7)
+
+ def _Loss7(m,i):
+ return m.loss *m.F[i]*0.435 *m.CODfeedactual
+ m.Loss7=Expression(m.t,rule=_Loss7)
+
+ def _DWN8(m,i):
+ return (1-m.loss)* m.F[i]* 0.204 *m.CODsolidactual
+ m.DWN8=Expression(m.t,rule=_DWN8)
+
+ def _Loss8(m,i):
+ return m.loss *m.F[i]*0.145 *m.CODfeedactual
+ m.Loss8=Expression(m.t,rule=_Loss8)
+
+ def _DWN9(m,i):
+ return (1-m.loss)* m.F[i]* 0.39 *m.CODsolidactual
+ m.DWN9=Expression(m.t,rule=_DWN9)
+
+ def _Loss9(m,i):
+ return m.loss *m.F[i]*0.275 *m.CODfeedactual
+ m.Loss9=Expression(m.t,rule=_Loss9)
+
+ def _r1(m,i):
+ return m.muh * (m.z6[i]/(m.KS+m.z6[i]))*(m.z3[i]/(m.KOH+m.z3[i])) * m.z10[i]
+ m.r1=Expression(m.t,rule=_r1)
+
+ def _r1b(m,i):
+ return m.muh * (m.z6[i]/(m.KS+m.z6[i]))*(m.z3[i]/(m.KOH+m.z3[i])) * m.z14[i]
+ m.r1b=Expression(m.t,rule=_r1b)
+
+
+ #On NO2
+ def _r2(m,i):
+ return m.muh*m.etano2*(m.z6[i]/(m.KS+m.z6[i])) * (m.KOH/(m.KOH + m.z3[i])) * (m.z61[i]/(m.KNOx + m.z61[i])) * (m.z61[i]/(m.z12[i] + m.z61[i])) * m.z10[i]
+ m.r2=Expression(m.t,rule=_r2)
+
+ def _r2b(m,i):
+ return m.muh*m.etano2*(m.z6[i]/(m.KS+m.z6[i])) * (m.KOH/(m.KOH + m.z3[i])) * (m.z61[i]/(m.KNOx + m.z61[i])) * (m.z61[i]/(m.z12[i] + m.z61[i]))* m.z14[i]
+ m.r2b=Expression(m.t,rule=_r2b)
+
+
+ #On NO3
+ def _rr2(m,i):
+ return m.muh*m.etano3*(m.z6[i]/(m.KS+m.z6[i])) * (m.KOH/(m.KOH + m.z3[i])) * (m.z12[i]/(m.KNOx + m.z12[i])) * (m.z12[i]/(m.z12[i] + m.z61[i])) * m.z10[i]
+ m.rr2=Expression(m.t,rule=_rr2)
+
+ def _rr2b(m,i):
+ return m.muh*m.etano3*(m.z6[i]/(m.KS+m.z6[i])) * (m.KOH/(m.KOH + m.z3[i])) * (m.z12[i]/(m.KNOx + m.z12[i])) * (m.z12[i]/(m.z12[i] + m.z61[i]))* m.z14[i]
+ m.rr2b=Expression(m.t,rule=_rr2b)
+
+
+ #growth of AOB
+ def _r3(m,i):
+ return m.muAOB*(m.z0[i]/(m.KNH+m.z0[i])) * (m.z3[i]/(m.KOA + m.z3[i])) * m.z11[i]
+ m.r3=Expression(m.t,rule=_r3)
+
+ def _r3b(m,i):
+ return m.muAOB*(m.z0[i]/(m.KNH+m.z0[i])) * (m.z3[i]/(m.KOA + m.z3[i])) * m.z15[i]
+ m.r3b=Expression(m.t,rule=_r3b)
+
+
+ #growth of NOB
+ def _rr3(m,i):
+ return m.muNOB*(m.z61[i]/(m.KNO2+m.z61[i])) * (m.z3[i]/(m.KOA + m.z3[i])) * (m.KNHI/(m.KNHI + m.z0[i]))* m.z64[i]
+ m.rr3=Expression(m.t,rule=_rr3)
+
+ def _rr3b(m,i):
+ return m.muNOB*(m.z61[i]/(m.KNO2+m.z61[i])) * (m.z3[i]/(m.KOA + m.z3[i])) * (m.KNHI/(m.KNHI + m.z0[i]))* m.z67[i]
+ m.rr3b=Expression(m.t,rule=_rr3b)
+
+
+ def _r4(m,i):
+ return m.bH*m.z10[i]
+ m.r4=Expression(m.t,rule=_r4)
+
+ def _r4b(m,i):
+ return m.bH*m.z14[i]
+ m.r4b=Expression(m.t,rule=_r4b)
+
+
+ def _r5(m,i):
+ return m.bAOB*m.z11[i]
+ m.r5=Expression(m.t,rule=_r5)
+
+ def _r5b(m,i):
+ return m.bAOB*m.z15[i]
+ m.r5b=Expression(m.t,rule=_r5b)
+
+ def _rr5(m,i):
+ return m.bNOB*m.z64[i]
+ m.rr5=Expression(m.t,rule=_rr5)
+
+ def _rr5b(m,i):
+ return m.bNOB*m.z67[i]
+ m.rr5b=Expression(m.t,rule=_rr5b)
+
+ def _r6(m,i):
+ return m.Ka* m.z1[i]* (m.z10[i] + m.z14[i])
+ m.r6=Expression(m.t,rule=_r6)
+
+ def _r7(m,i):
+ return m.kh*((m.z7[i]/m.z10[i])/(m.Kx+(m.z7[i]/m.z10[i])) ) * ( (m.z3[i]/(m.KOH+m.z3[i])) + m.etag*(m.KOH/ (m.KOH + m.z3[i]))* ((m.z12[i]+m.z61[i])/(m.KNOx+m.z12[i]+m.z61[i]) ) )*m.z10[i]
+ m.r7=Expression(m.t,rule=_r7)
+
+ def _r7b(m,i):
+ return m.kh*((m.z16[i]/m.z14[i])/(m.Kx+(m.z16[i]/m.z14[i])) ) * ( (m.z3[i]/(m.KOH+m.z3[i])) + m.etag*(m.KOH/ (m.KOH + m.z3[i]))* ((m.z12[i]+m.z61[i])/(m.KNOx+m.z12[i]+m.z61[i]) ) )*m.z14[i]
+ m.r7b=Expression(m.t,rule=_r7b)
+
+ def _r8(m,i):
+ return m.r7[i]*(m.z2[i]/m.z7[i])
+ m.r8=Expression(m.t,rule=_r8)
+
+ def _r8b(m,i):
+ return m.r7b[i]*(m.z17[i]/m.z16[i])
+ m.r8b=Expression(m.t,rule=_r8b)
+
+ def _r9(m,i):
+ return m.kka*((m.z9[i]+m.z7[i]+m.z10[i]+m.z11[i]+m.z64[i]+m.z13[i])/m.icv+m.a)*m.z9[i]
+ m.r9=Expression(m.t,rule=_r9)
+
+ def _r10(m,i):
+ return m.kka*((m.z9[i]+m.z7[i]+m.z10[i]+m.z11[i]+m.z64[i]+m.z13[i])/m.icv+m.a)*m.z7[i]
+ m.r10=Expression(m.t,rule=_r10)
+
+ def _r11(m,i):
+ return m.kka*((m.z9[i]+m.z7[i]+m.z10[i]+m.z11[i]+m.z64[i]+m.z13[i])/m.icv+m.a)*m.z10[i]
+ m.r11=Expression(m.t,rule=_r11)
+
+ def _r12(m,i):
+ return m.kka*((m.z9[i]+m.z7[i]+m.z10[i]+m.z11[i]+m.z64[i]+m.z13[i])/m.icv+m.a)*m.z11[i]
+ m.r12=Expression(m.t,rule=_r12)
+
+ def _rr12(m,i):
+ return m.kka*((m.z9[i]+m.z7[i]+m.z10[i]+m.z11[i]+m.z64[i]+m.z13[i])/m.icv+m.a)*m.z64[i]
+ m.rr12=Expression(m.t,rule=_rr12)
+
+ def _r13(m,i):
+ return m.kka*((m.z9[i]+m.z7[i]+m.z10[i]+m.z11[i]+m.z64[i]+m.z13[i])/m.icv+m.a)*m.z13[i]
+ m.r13=Expression(m.t,rule=_r13)
+
+ def _r14(m,i):
+ return m.kka*((m.z9[i]+m.z7[i]+m.z10[i]+m.z11[i]+m.z64[i]+m.z13[i])/m.icv+m.a)*m.z2[i]
+ m.r14=Expression(m.t,rule=_r14)
+
+ def _r15(m,i):
+ return m.kkd*((m.z18[i]+m.z16[i]+m.z14[i]+m.z15[i]+m.z67[i]+m.z19[i])/m.icv+m.a)*m.z18[i]
+ m.r15=Expression(m.t,rule=_r15)
+
+ def _r16(m,i):
+ return m.kkd*((m.z18[i]+m.z16[i]+m.z14[i]+m.z15[i]+m.z67[i]+m.z19[i])/m.icv+m.a)*m.z16[i]
+ m.r16=Expression(m.t,rule=_r16)
+
+ def _r17(m,i):
+ return m.kkd*((m.z18[i]+m.z16[i]+m.z14[i]+m.z15[i]+m.z67[i]+m.z19[i])/m.icv+m.a)*m.z14[i]
+ m.r17=Expression(m.t,rule=_r17)
+
+ def _r18(m,i):
+ return m.kkd*((m.z18[i]+m.z16[i]+m.z14[i]+m.z15[i]+m.z67[i]+m.z19[i])/m.icv+m.a)*m.z15[i]
+ m.r18=Expression(m.t,rule=_r18)
+
+ def _rr18(m,i):
+ return m.kkd*((m.z18[i]+m.z16[i]+m.z14[i]+m.z15[i]+m.z67[i]+m.z19[i])/m.icv+m.a)*m.z67[i]
+ m.rr18=Expression(m.t,rule=_rr18)
+
+ def _r19(m,i):
+ return m.kkd*((m.z18[i]+m.z16[i]+m.z14[i]+m.z15[i]+m.z67[i]+m.z19[i])/m.icv+m.a)*m.z19[i]
+ m.r19=Expression(m.t,rule=_r19)
+
+ def _r20(m,i):
+ return m.kkd*((m.z18[i]+m.z16[i]+m.z14[i]+m.z15[i]+m.z67[i]+m.z19[i])/m.icv+m.a)*m.z17[i]
+ m.r20=Expression(m.t,rule=_r20)
+
+
+ def _bgro(m,i):
+ sk1=m.z69[i]*0.006
+ sk2=log(sk1)
+ return -0.7*sk2
+ m.bgro=Expression(m.t,rule=_bgro)
+
+
+ ############ reaction terms ############
+
+ #Readily biodegradable organic matter
+ def _rs(m,i):
+ return -(1/m.YH)*(m.r1[i]+m.r2[i]+m.r1b[i]+m.r2b[i]+m.rr2[i]+m.rr2b[i]) + m.r7[i]+m.r7b[i]
+ m.rs=Expression(m.t,rule=_rs)
+
+
+ #Inert organic matter particulate
+ def _rIX(m,i):
+ return -m.r9[i] + m.r15[i]
+ m.rIX=Expression(m.t,rule=_rIX)
+
+
+ #Slowly biodegradable organic matter
+ def _rX(m,i):
+ return (1-m.fp)*(m.r4[i]+m.r5[i]+m.rr5[i]) -m.r7[i] -m.r10[i]+ m.r16[i]
+ m.rX=Expression(m.t,rule=_rX)
+
+
+ #Heterotrphps
+ def _rBH(m,i):
+ return m.r1[i]+m.r2[i]+m.rr2[i]-m.r4[i] -m.r11[i] + m.r17[i]
+ m.rBH=Expression(m.t,rule=_rBH)
+
+
+
+ #AOB
+ def _rAOB(m,i):
+ return m.r3[i]-m.r5[i]-m.r12[i]+m.r18[i]
+ m.rAOB=Expression(m.t,rule=_rAOB)
+
+
+ #NOB
+ def _rNOB(m,i):
+ return m.rr3[i]-m.rr5[i]-m.rr12[i]+m.rr18[i]
+ m.rNOB=Expression(m.t,rule=_rNOB)
+
+
+ #Particulate products from biomass deccay
+ def _rP(m,i):
+ return m.fp*(m.r4[i]+m.r5[i]+m.rr5[i]) - m.r13[i] + m.r19[i]
+ m.rP=Expression(m.t,rule=_rP)
+
+
+ #Inert particulate biofilm
+ def _rIXb(m,i):
+ return m.r9[i]-m.r15[i]
+ m.rIXb=Expression(m.t,rule=_rIXb)
+
+
+ #Slowly biodegradable organic matter biofilm
+ def _rXb(m,i):
+ return (1-m.fp)*(m.r4b[i]+m.r5b[i]+m.rr5b[i]) -m.r7b[i] +m.r10[i] -m.r16[i]
+ m.rXb=Expression(m.t,rule=_rXb)
+
+
+ #Heterotrphps biofilm
+ def _rBHb(m,i):
+ return m.r1b[i]+m.r2b[i]+m.rr2b[i]-m.r4b[i] +m.r11[i] - m.r17[i]
+ m.rBHb=Expression(m.t,rule=_rBHb)
+
+
+ #AOB biofilm
+ def _rAOBb(m,i):
+ return m.r3b[i]-m.r5b[i]+m.r12[i]-m.r18[i]
+ m.rAOBb=Expression(m.t,rule=_rAOBb)
+
+ #NOB
+ def _rNOBb(m,i):
+ return m.rr3b[i]-m.rr5b[i]+m.rr12[i]-m.rr18[i]
+ m.rNOBb=Expression(m.t,rule=_rNOBb)
+
+
+ #Particulate products from biomass deccay biofilm
+ def _rPb(m,i):
+ return m.fp*(m.r4b[i]+m.r5b[i]+m.rr5b[i]) + m.r13[i] - m.r19[i]
+ m.rPb=Expression(m.t,rule=_rPb)
+
+
+ #Nitrite
+ def _rNO2(m,i):
+ return ((1-1/m.YH)/(1.72))*(m.r2[i] + m.r2b[i]) + (1/m.YAOB)*(m.r3[i] +m.r3b[i])- (1/m.YNOB)*(m.rr3[i] +m.rr3b[i])
+ m.rNO2=Expression(m.t,rule=_rNO2)
+
+
+ #Nitrate
+ def _rNO3(m,i):
+ return ((1-1/m.YH)/(2.86))*(m.rr2[i] + m.rr2b[i]) + (1/m.YNOB)*(m.rr3[i] +m.rr3b[i])
+ m.rNO3=Expression(m.t,rule=_rNO3)
+
+
+ # Ammonia
+ def _ramm(m,i):
+ return -m.iXB*(m.r1[i]+m.r1b[i]+ m.r2[i] + m.r2b[i]+m.rr2[i] + m.rr2b[i]) - (m.iXB + (1/m.YAOB)) * (m.r3[i] + m.r3b[i]) - m.iXB*(m.rr3[i] + m.rr3b[i]) + m.r6[i]
+ m.ramm=Expression(m.t,rule=_ramm)
+
+
+
+ #soluble organic nitrogen
+ def _rSND(m,i):
+ return -m.r6[i] + m.r8[i] + m.r8b[i]
+ m.rSND=Expression(m.t,rule=_rSND)
+
+
+
+ #particulate organic nitrogen
+ def _rXND(m,i):
+ return (m.iXB-m.fp*m.iXP)*(m.r4[i]+m.r5[i]+m.rr5[i]) -m.r8[i] -m.r14[i] + m.r20[i]
+ m.rXND=Expression(m.t,rule=_rXND)
+
+
+
+ #particulate organic nitrogen biofilm
+ def _rXNDb(m,i):
+ return (m.iXB-m.fp*m.iXP)*(m.r4b[i]+m.r5b[i]+m.rr5b[i]) -m.r8b[i] + m.r14[i] - m.r20[i]
+ m.rXNDb=Expression(m.t,rule=_rXNDb)
+
+
+ ###### second reactor
+
+ # Reactions in the first reactor
+ def _rec1(m,i):
+ return m.muh * (m.z26[i]/(m.KS+m.z26[i]))*(m.z23[i]/(m.KOH+m.z23[i])) * m.z30[i]
+ m.rec1=Expression(m.t,rule=_rec1)
+
+ def _rec1b(m,i):
+ return m.muh * (m.z26[i]/(m.KS+m.z26[i]))*(m.z23[i]/(m.KOH+m.z23[i])) * m.z34[i]
+ m.rec1b=Expression(m.t,rule=_rec1b)
+
+
+ #On NO2
+ def _rec2(m,i):
+ return m.muh*m.etano2*(m.z26[i]/(m.KS+m.z26[i])) * (m.KOH/(m.KOH + m.z23[i])) * (m.z62[i]/(m.KNOx + m.z62[i])) * (m.z62[i]/(m.z32[i] + m.z62[i]))* m.z30[i]
+ m.rec2=Expression(m.t,rule=_rec2)
+
+ def _rec2b(m,i):
+ return m.muh*m.etano2*(m.z26[i]/(m.KS+m.z26[i])) * (m.KOH/(m.KOH + m.z23[i])) * (m.z62[i]/(m.KNOx + m.z62[i])) * (m.z62[i]/(m.z32[i] + m.z62[i])) * m.z34[i]
+
+ m.rec2b=Expression(m.t,rule=_rec2b)
+
+
+ #On NO3
+ def _recrec2(m,i):
+ return m.muh*m.etano3*(m.z26[i]/(m.KS+m.z26[i])) * (m.KOH/(m.KOH + m.z23[i])) * (m.z32[i]/(m.KNOx + m.z32[i])) * (m.z32[i]/(m.z32[i] + m.z62[i]))* m.z30[i]
+ m.recrec2=Expression(m.t,rule=_recrec2)
+
+ def _recrec2b(m,i):
+ return m.muh*m.etano3*(m.z26[i]/(m.KS+m.z26[i])) * (m.KOH/(m.KOH + m.z23[i])) * (m.z32[i]/(m.KNOx + m.z32[i])) * (m.z32[i]/(m.z32[i] + m.z62[i])) * m.z34[i]
+ m.recrec2b=Expression(m.t,rule=_recrec2b)
+
+
+ #growth of AOB
+ def _rec3(m,i):
+ return m.muAOB*(m.z20[i]/(m.KNH+m.z20[i])) * (m.z23[i]/(m.KOA + m.z23[i])) * m.z31[i]
+ m.rec3=Expression(m.t,rule=_rec3)
+
+ def _rec3b(m,i):
+ return m.muAOB*(m.z20[i]/(m.KNH+m.z20[i])) * (m.z23[i]/(m.KOA + m.z23[i])) * m.z35[i]
+ m.rec3b=Expression(m.t,rule=_rec3b)
+
+
+ #growth of NOB
+
+ def _recrec3(m,i):
+ return m.muNOB*(m.z62[i]/(m.KNO2+m.z62[i])) * (m.z23[i]/(m.KOA + m.z23[i])) * (m.KNHI/(m.KNHI + m.z20[i]))* m.z65[i]
+ m.recrec3=Expression(m.t,rule=_recrec3)
+
+ def _recrec3b(m,i):
+ return m.muNOB*(m.z62[i]/(m.KNO2+m.z62[i])) * (m.z23[i]/(m.KOA + m.z23[i])) * (m.KNHI/(m.KNHI + m.z20[i]))* m.z68[i]
+ m.recrec3b=Expression(m.t,rule=_recrec3b)
+
+
+ def _rec4(m,i):
+ return m.bH*m.z30[i]
+ m.rec4=Expression(m.t,rule=_rec4)
+
+ def _rec4b(m,i):
+ return m.bH*m.z34[i]
+ m.rec4b=Expression(m.t,rule=_rec4b)
+
+
+ def _rec5(m,i):
+ return m.bAOB*m.z31[i]
+ m.rec5=Expression(m.t,rule=_rec5)
+
+ def _rec5b(m,i):
+ return m.bAOB*m.z35[i]
+ m.rec5b=Expression(m.t,rule=_rec5b)
+
+ def _recrec5(m,i):
+ return m.bNOB*m.z65[i]
+ m.recrec5=Expression(m.t,rule=_recrec5)
+
+ def _recrec5b(m,i):
+ return m.bNOB*m.z68[i]
+ m.recrec5b=Expression(m.t,rule=_recrec5b)
+
+ def _rec6(m,i):
+ return m.Ka* m.z21[i]* (m.z30[i] + m.z34[i])
+ m.rec6=Expression(m.t,rule=_rec6)
+
+ def _rec7(m,i):
+ return m.kh*((m.z27[i]/m.z30[i])/(m.Kx+(m.z27[i]/m.z30[i])) ) * ( (m.z23[i]/(m.KOH+m.z23[i])) + m.etag*(m.KOH/ (m.KOH + m.z23[i]))* ((m.z32[i]+m.z62[i])/(m.KNOx+m.z32[i]+m.z62[i]) ) )*m.z30[i]
+ m.rec7=Expression(m.t,rule=_rec7)
+
+ def _rec7b(m,i):
+ return m.kh*((m.z36[i]/m.z34[i])/(m.Kx+(m.z36[i]/m.z34[i])) ) * ( (m.z23[i]/(m.KOH+m.z23[i])) + m.etag*(m.KOH/ (m.KOH + m.z23[i]))* ((m.z32[i]+m.z62[i])/(m.KNOx+m.z32[i]+m.z62[i]) ) )*m.z34[i]
+ m.rec7b=Expression(m.t,rule=_rec7b)
+
+ def _rec8(m,i):
+ return m.rec7[i]*(m.z22[i]/m.z27[i])
+ m.rec8=Expression(m.t,rule=_rec8)
+
+ def _rec8b(m,i):
+ return m.rec7b[i]*(m.z37[i]/m.z36[i])
+ m.rec8b=Expression(m.t,rule=_rec8b)
+
+ def _rec9(m,i):
+ return m.kka*((m.z29[i]+m.z27[i]+m.z30[i]+m.z31[i]+m.z65[i]+m.z33[i])/m.icv+m.a)*m.z29[i]
+ m.rec9=Expression(m.t,rule=_rec9)
+
+ def _rec10(m,i):
+ return m.kka*((m.z29[i]+m.z27[i]+m.z30[i]+m.z31[i]+m.z65[i]+m.z33[i])/m.icv+m.a)*m.z27[i]
+ m.rec10=Expression(m.t,rule=_rec10)
+
+ def _rec11(m,i):
+ return m.kka*((m.z29[i]+m.z27[i]+m.z30[i]+m.z31[i]+m.z65[i]+m.z33[i])/m.icv+m.a)*m.z30[i]
+ m.rec11=Expression(m.t,rule=_rec11)
+
+ def _rec12(m,i):
+ return m.kka*((m.z29[i]+m.z27[i]+m.z30[i]+m.z31[i]+m.z65[i]+m.z33[i])/m.icv+m.a)*m.z31[i]
+ m.rec12=Expression(m.t,rule=_rec12)
+
+ def _recrec12(m,i):
+ return m.kka*((m.z29[i]+m.z27[i]+m.z30[i]+m.z31[i]+m.z65[i]+m.z33[i])/m.icv+m.a)*m.z65[i]
+ m.recrec12=Expression(m.t,rule=_recrec12)
+
+ def _rec13(m,i):
+ return m.kka*((m.z29[i]+m.z27[i]+m.z30[i]+m.z31[i]+m.z65[i]+m.z33[i])/m.icv+m.a)*m.z33[i]
+ m.rec13=Expression(m.t,rule=_rec13)
+
+ def _rec14(m,i):
+ return m.kka*((m.z29[i]+m.z27[i]+m.z30[i]+m.z31[i]+m.z65[i]+m.z33[i])/m.icv+m.a)*m.z22[i]
+ m.rec14=Expression(m.t,rule=_rec14)
+
+ def _rec15(m,i):
+ return m.kkd*((m.z38[i]+m.z36[i]+m.z34[i]+m.z35[i]+m.z68[i]+m.z39[i])/m.icv+m.a)*m.z38[i]
+ m.rec15=Expression(m.t,rule=_rec15)
+
+ def _rec16(m,i):
+ return m.kkd*((m.z38[i]+m.z36[i]+m.z34[i]+m.z35[i]+m.z68[i]+m.z39[i])/m.icv+m.a)*m.z36[i]
+ m.rec16=Expression(m.t,rule=_rec16)
+
+ def _rec17(m,i):
+ return m.kkd*((m.z38[i]+m.z36[i]+m.z34[i]+m.z35[i]+m.z68[i]+m.z39[i])/m.icv+m.a)*m.z34[i]
+ m.rec17=Expression(m.t,rule=_rec17)
+
+ def _rec18(m,i):
+ return m.kkd*((m.z38[i]+m.z36[i]+m.z34[i]+m.z35[i]+m.z68[i]+m.z39[i])/m.icv+m.a)*m.z35[i]
+ m.rec18=Expression(m.t,rule=_rec18)
+
+ def _recrec18(m,i):
+ return m.kkd*((m.z38[i]+m.z36[i]+m.z34[i]+m.z35[i]+m.z68[i]+m.z39[i])/m.icv+m.a)*m.z68[i]
+ m.recrec18=Expression(m.t,rule=_recrec18)
+
+ def _rec19(m,i):
+ return m.kkd*((m.z38[i]+m.z36[i]+m.z34[i]+m.z35[i]+m.z68[i]+m.z39[i])/m.icv+m.a)*m.z39[i]
+ m.rec19=Expression(m.t,rule=_rec19)
+
+ def _rec20(m,i):
+ return m.kkd*((m.z38[i]+m.z36[i]+m.z34[i]+m.z35[i]+m.z68[i]+m.z39[i])/m.icv+m.a)*m.z37[i]
+ m.rec20=Expression(m.t,rule=_rec20)
+
+
+
+ ############ reaction terms ############
+
+ #Readily biodegradable organic matter
+ def _recs(m,i):
+ return -(1/m.YH)*(m.rec1[i]+m.rec2[i]+m.rec1b[i]+m.rec2b[i]+m.recrec2[i]+m.recrec2b[i]) + m.rec7[i]+m.rec7b[i]
+ m.recs=Expression(m.t,rule=_recs)
+
+
+
+ #Inert organic matter particulate
+ def _recIX(m,i):
+ return -m.rec9[i] + m.rec15[i]
+ m.recIX=Expression(m.t,rule=_recIX)
+
+
+ #Slowly biodegradable organic matter
+ def _recX(m,i):
+ return (1-m.fp)*(m.rec4[i]+m.rec5[i]+m.recrec5[i]) -m.rec7[i] -m.rec10[i]+ m.rec16[i]
+
+ m.recX=Expression(m.t,rule=_recX)
+
+
+
+ #Heterotrphps
+ def _recBH(m,i):
+ return m.rec1[i]+m.rec2[i]+m.recrec2[i]-m.rec4[i] -m.rec11[i] + m.rec17[i]
+ m.recBH=Expression(m.t,rule=_recBH)
+
+
+ #AOB
+ def _recAOB(m,i):
+ return m.rec3[i]-m.rec5[i]-m.rec12[i]+m.rec18[i]
+ m.recAOB=Expression(m.t,rule=_recAOB)
+
+
+
+ #NOB
+ def _recNOB(m,i):
+ return m.recrec3[i]-m.recrec5[i]-m.recrec12[i]+m.recrec18[i]
+ m.recNOB=Expression(m.t,rule=_recNOB)
+
+
+ #Particulate products from biomass deccay
+ def _recP(m,i):
+ return m.fp*(m.rec4[i]+m.rec5[i]+m.recrec5[i]) - m.rec13[i] + m.rec19[i]
+ m.recP=Expression(m.t,rule=_recP)
+
+
+ #Inert particulate biofilm
+ def _recIXb(m,i):
+ return m.rec9[i]-m.rec15[i]
+ m.recIXb=Expression(m.t,rule=_recIXb)
+
+
+
+ #Slowly biodegradable organic matter biofilm
+ def _recXb(m,i):
+ return (1-m.fp)*(m.rec4b[i]+m.rec5b[i]+m.recrec5b[i]) -m.rec7b[i] +m.rec10[i] -m.rec16[i]
+ m.recXb=Expression(m.t,rule=_recXb)
+
+
+ #Heterotrphps biofilm
+ def _recBHb(m,i):
+ return m.rec1b[i]+m.rec2b[i]+m.recrec2b[i]-m.rec4b[i] +m.rec11[i] - m.rec17[i]
+ m.recBHb=Expression(m.t,rule=_recBHb)
+
+
+ #AOB biofilm
+ def _recAOBb(m,i):
+ return m.rec3b[i]-m.rec5b[i]+m.rec12[i]-m.rec18[i]
+ m.recAOBb=Expression(m.t,rule=_recAOBb)
+
+
+ #NOB
+ def _recNOBb(m,i):
+ return m.recrec3b[i]-m.recrec5b[i]+ m.recrec12[i]- m.recrec18[i]
+ m.recNOBb=Expression(m.t,rule=_recNOBb)
+
+
+ #Particulate products from biomass deccay biofilm
+ def _recPb(m,i):
+ return m.fp*(m.rec4b[i]+m.rec5b[i]+m.recrec5b[i]) + m.rec13[i] - m.rec19[i]
+ m.recPb=Expression(m.t,rule=_recPb)
+
+
+ def _rO(m,i):
+ return (1-(1/m.YH))*(m.r1[i] + m.r1b[i]) - ( (3.43-m.YAOB)/m.YAOB)*(m.r3[i] + m.r3b[i]) - ( (1.14-m.YNOB)/m.YNOB)*(m.rr3[i] + m.rr3b[i])
+ m.rO=Expression(m.t,rule=_rO)
+
+ def _recO(m,i):
+ return (1-(1/m.YH))*(m.rec1[i] + m.rec1b[i]) - ( (3.43-m.YAOB)/m.YAOB)*(m.rec3[i] + m.rec3b[i]) - ( (1.14-m.YNOB)/m.YNOB)*(m.recrec3[i] +m.recrec3b[i])
+ m.recO=Expression(m.t,rule=_recO)
+
+
+ #Nitrite
+ def _recNO2(m,i):
+ return ((1-1/m.YH)/(1.72))*(m.rec2[i] + m.rec2b[i]) + (1/m.YAOB)*(m.rec3[i] +m.rec3b[i])- (1/m.YNOB)*(m.recrec3[i] +m.recrec3b[i])
+ m.recNO2=Expression(m.t,rule=_recNO2)
+
+
+ #Nitrate
+ def _recNO3(m,i):
+ return ((1-1/m.YH)/(2.86))*(m.recrec2[i] + m.recrec2b[i]) + (1/m.YNOB)*(m.recrec3[i] +m.recrec3b[i])
+ m.recNO3=Expression(m.t,rule=_recNO3)
+
+
+ # Ammonia
+ def _recamm(m,i):
+ return -m.iXB*(m.rec1[i]+m.rec1b[i]+ m.rec2[i] + m.rec2b[i]+m.recrec2[i] + m.recrec2b[i]) - (m.iXB + (1/m.YAOB)) * (m.rec3[i] + m.rec3b[i]) - m.iXB*(m.recrec3[i] + m.recrec3b[i]) + m.rec6[i]
+ m.recamm=Expression(m.t,rule=_recamm)
+
+
+ #soluble organic nitrogen
+ def _recSND(m,i):
+ return -m.rec6[i] + m.rec8[i] + m.rec8b[i]
+ m.recSND=Expression(m.t,rule=_recSND)
+
+
+ #particulate organic nitrogen
+ def _recXND(m,i):
+ return (m.iXB-m.fp*m.iXP)*(m.rec4[i]+m.rec5[i]+m.recrec5[i]) -m.rec8[i] -m.rec14[i] + m.rec20[i]
+ m.recXND=Expression(m.t,rule=_recXND)
+
+
+ #particulate organic nitrogen biofilm
+ def _recXNDb(m,i):
+ return (m.iXB-m.fp*m.iXP)*(m.rec4b[i]+m.rec5b[i]+m.recrec5b[i]) -m.rec8b[i] + m.rec14[i] - m.rec20[i]
+ m.recXNDb=Expression(m.t,rule=_recXNDb)
+
+
+ #### fish tanks ########
+ def _c1(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni54=1+value(m.ni54)
+
+ return m.dx1dt[i] +(1/m.tav1)*m.z54[i] - (1/m.tav1)*m.DWN1[i] ==0
+ m.c1=Constraint(m.t, rule=_c1)
+
+
+ def _c2(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni40=1+value(m.ni40)
+ return m.dy1dt[i] - ((m.Q1)/m.V)*(0.5*m.z20[i]+0.5*m.z0[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z40[i]) - (1/m.V)* (m.z54[i]) ==0
+ m.c2=Constraint(m.t, rule=_c2)
+
+ def _c3(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni55=1+value(m.ni55)
+ return m.dx2dt[i] +(1/m.tav1)*m.z55[i]- (1/m.tav1)*m.DWN2[i] ==0
+ m.c3=Constraint(m.t, rule=_c3)
+
+
+ def _c4(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni41=1+value(m.ni41)
+ return m.dy2dt[i] - ((m.Q1)/m.V)*(0.5*m.z21[i]+0.5*m.z1[i]) +((m.Q1+m.Q4[i])/m.V)*(m.z41[i]) - (1/m.V)* (m.z55[i]) - (1/m.V)* (m.Loss2[i] ) ==0
+ m.c4=Constraint(m.t, rule=_c4)
+
+
+ def _c5(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni56=1+value(m.ni56)
+ return m.dx3dt[i] +(1/m.tav1)*m.z56[i]- (1/m.tav1)*m.DWN3[i]==0
+ m.c5=Constraint(m.t, rule=_c5)
+
+ def _c6(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni42=1+value(m.ni42)
+ return m.dy3dt[i] - ((m.Q1)/m.V)*(0.5*m.z22[i]+0.5*m.z2[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z42[i]) - (1/m.V)* (m.z56[i]) - (1/m.V)* (m.Loss3[i] ) ==0
+ m.c6=Constraint(m.t, rule=_c6)
+
+
+ def _c7(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni57=1+value(m.ni57)
+ return m.dx6dt[i] + (1/m.tav1)*m.z57[i]- (1/m.tav1)*m.DWN6[i] ==0
+ m.c7=Constraint(m.t, rule=_c7)
+
+
+ def _c8(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni45=1+value(m.ni45)
+ return m.dy6dt[i] - ((m.Q1)/m.V)*(0.5*m.z26[i]+0.5*m.z6[i]) +((m.Q1+m.Q4[i])/m.V)*(m.z45[i]) - (1/m.V)* (m.z57[i]) - (1/m.V)* (m.Loss6[i] ) ==0
+ m.c8=Constraint(m.t, rule=_c8)
+
+
+ def _c9(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni58=1+value(m.ni58)
+ return m.dx7dt[i] +(1/m.tav1)*m.z58[i]- (1/m.tav1)*m.DWN7[i] ==0
+ m.c9=Constraint(m.t, rule=_c9)
+
+
+ def _c10(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni46=1+value(m.ni46)
+ return m.dy7dt[i] - ((m.Q1)/m.V)*(0.5*m.z27[i]+0.5*m.z7[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z46[i]) - (1/m.V)* (m.z58[i]) - (1/m.V)* (m.Loss7[i] ) ==0
+ m.c10=Constraint(m.t, rule=_c10)
+
+ def _c11(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni59=1+value(m.ni59)
+ return m.dx8dt[i] + (1/m.tav1)*m.z59[i] - (1/m.tav1)*m.DWN8[i] ==0
+ m.c11=Constraint(m.t, rule=_c11)
+
+
+ def _c12(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni47=1+value(m.ni47)
+# m.ni47+=1
+ return m.dy8dt[i] - ((m.Q1)/m.V)*(0.5*m.z28[i]+0.5*m.z8[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z47[i]) - (1/m.V)* (m.z59[i]) - (1/m.V)* (m.Loss8[i] )==0
+ m.c12=Constraint(m.t, rule=_c12)
+
+ def _c13(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+# m.ni60+=1
+ m.ni60=1+value(m.ni60)
+ return m.dx9dt[i] + (1/m.tav1)*m.z60[i] - (1/m.tav1)*m.DWN9[i] ==0
+ m.c13=Constraint(m.t, rule=_c13)
+
+ def _c14(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni48=1+value(m.ni48)
+ return m.dy9dt[i] - ((m.Q1)/m.V)*(0.5*m.z29[i]+0.5*m.z9[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z48[i]) - (1/m.V)* (m.z60[i]) - (1/m.V)* (m.Loss9[i] ) ==0
+
+ m.c14=Constraint(m.t, rule=_c14)
+
+
+
+ def _c15(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni50=1+value(m.ni50)
+ return m.dy11dt[i] - ((m.Q1)/m.V)*(0.5*m.z33[i]+0.5*m.z13[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z50[i]) ==0
+ m.c15=Constraint(m.t, rule=_c15)
+
+ #AOB
+ def _c16(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni51=1+value(m.ni51)
+ return m.dy12dt[i] - ((m.Q1)/m.V)*(0.5*m.z31[i]+0.5*m.z11[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z51[i]) - ((m.Q)/m.V)*(m.Aint) ==0
+
+ m.c16=Constraint(m.t, rule=_c16)
+
+
+ #NOB
+
+ def _c17(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni66=1+value(m.ni66)
+ return m.ddy12dt[i] - ((m.Q1)/m.V)*(0.5*m.z65[i]+0.5*m.z64[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z66[i]) - ((m.Q)/m.V)*(m.Aint) ==0
+ m.c17=Constraint(m.t, rule=_c17)
+
+ def _c18(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni52=1+value(m.ni52)
+ return m.dy13dt[i] - ((m.Q1)/m.V)*(0.5*m.z30[i]+0.5*m.z10[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z52[i]) - ((m.Q)/m.V)*(m.Hint) ==0
+ m.c18=Constraint(m.t, rule=_c18)
+
+ #NO3
+ def _c19(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni53=1+value(m.ni53)
+ return m.dy14dt[i] - ((m.Q1)/m.V)*(0.5*m.z32[i]+0.5*m.z12[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z53[i]) ==0
+ m.c19=Constraint(m.t, rule=_c19)
+
+
+ #NO2
+ def _c20(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni63=1+value(m.ni63)
+ return m.ddy14dt[i] - ((m.Q1)/m.V)*(0.5*m.z62[i]+0.5*m.z61[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z63[i]) ==0
+ m.c20=Constraint(m.t, rule=_c20)
+
+
+
+################# wastewater treatment ###################
+
+ #Inert soluble organic matter
+ def _cz1(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni63=1+value(m.ni63)
+ return m.dzSIdt[i] - (m.Q/m.Vw1)*(m.z47[i]-m.z8[i]) ==0
+ m.cz1=Constraint(m.t, rule=_cz1)
+
+
+ # readily organic soluble
+ def _cz2(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni6=1+value(m.ni6)
+ return m.dzSdt[i] - (m.Q/m.Vw1)*(m.z45[i]-m.z6[i]) - m.rs[i] ==0
+ m.cz2=Constraint(m.t, rule=_cz2)
+
+
+ #Inert particulate organic matter
+ def _cz3(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni9=1+value(m.ni9)
+ return m.dzXIdt[i] - (m.Q/m.Vw1)*(0.52*((m.Q1+m.Q4[i])/m.Q1)*m.z48[i]-m.z9[i]) - m.rIX[i] ==0
+ m.cz3=Constraint(m.t, rule=_cz3)
+
+
+ #slowly organic matter
+ def _cz4(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni7=1+value(m.ni7)
+ return m.dzXdt[i] - (m.Q/m.Vw1)*(0.52*((m.Q1+m.Q4[i])/m.Q1)*m.z46[i]-m.z7[i]) - m.rX[i] ==0
+ m.cz4=Constraint(m.t, rule=_cz4)
+
+
+ #Heterotrophic
+ def _cz5(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni10=1+value(m.ni10)
+ return m.dXHdt[i] - (m.Q/m.Vw1)*(m.z52[i]-m.z10[i]) - m.rBH[i] ==0
+ m.cz5=Constraint(m.t, rule=_cz5)
+
+ #AOB
+ def _cz6(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni11=1+value(m.ni11)
+ return m.dAdt[i] - (m.Q/m.Vw1)*(m.z51[i]-m.z11[i]) - m.rAOB[i] ==0
+ m.cz6=Constraint(m.t, rule=_cz6)
+
+
+ #NOB
+ def _cz7(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni64=1+value(m.ni64)
+ return m.ddAdt[i] - (m.Q/m.Vw1)*(m.z66[i]-m.z64[i]) - m.rNOB[i] ==0
+ m.cz7=Constraint(m.t, rule=_cz7)
+
+ #Biomass decay
+ def _cz8(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni13=1+value(m.ni13)
+ return m.dXPdt[i] - (m.Q/m.Vw1)*(0.52*((m.Q1+m.Q4[i])/m.Q1)*m.z50[i]-m.z13[i]) - m.rP[i] ==0
+ m.cz8=Constraint(m.t, rule=_cz8)
+
+
+ #Inert particulate organic matter biofilm
+ def _cz9(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni18=1+value(m.ni18)
+ return m.dzXIbdt[i] - m.rIXb[i] ==0
+ m.cz9=Constraint(m.t, rule=_cz9)
+
+
+ #slowly organic matter biofilm
+ def _cz10(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni16=1+value(m.ni16)
+ return m.dzXbdt[i] - m.rXb[i] ==0
+ m.cz10=Constraint(m.t, rule=_cz10)
+
+
+ #Heterotrophic biofilm
+ def _cz11(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni14=1+value(m.ni14)
+ return m.dXHbdt[i] - m.rBHb[i] ==0
+ m.cz11=Constraint(m.t, rule=_cz11)
+
+
+
+ #AOB biofilm
+ def _cz12(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni15=1+value(m.ni15)
+ return m.dAbdt[i] - m.rAOBb[i] ==0
+ m.cz12=Constraint(m.t, rule=_cz12)
+
+
+ #NOB biofilm
+ def _cz13(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni67=1+value(m.ni67)
+ return m.ddAbdt[i] - m.rNOBb[i] ==0
+ m.cz13=Constraint(m.t, rule=_cz13)
+
+
+
+ #Biomass decay biofilm
+ def _cz14(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni19=1+value(m.ni19)
+ return m.dXPbdt[i] - m.rPb[i] ==0
+ m.cz14=Constraint(m.t, rule=_cz14)
+
+
+
+ #NO3
+ def _cz16(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni12=1+value(m.ni12)
+ return m.dzNOdt[i] - (m.Q/m.Vw1)*(m.z53[i]-m.z12[i]) - m.rNO3[i] ==0
+ m.cz16=Constraint(m.t, rule=_cz16)
+
+ #NO2
+ def _cz17(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni61=1+value(m.ni61)
+ return m.ddzNOdt[i] - (m.Q/m.Vw1)*(m.z63[i]-m.z61[i]) - m.rNO2[i] ==0
+ m.cz17=Constraint(m.t, rule=_cz17)
+
+
+ #ammonia
+ def _cz18(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni0=1+value(m.ni0)
+ return m.dNHdt[i] - (m.Q/m.Vw1)* (m.z40[i] -m.z0[i]) - m.ramm[i] ==0
+ m.cz18=Constraint(m.t, rule=_cz18)
+
+
+
+ #organic nitrogen soluble
+ def _cz19(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni1=1+value(m.ni1)
+ return m.dSNDdt[i] - (m.Q/m.Vw1)* (m.z41[i] -m.z1[i]) - m.rSND[i] ==0
+ m.cz19=Constraint(m.t, rule=_cz19)
+
+
+
+ #organic nitrogen particulate
+ def _cz20(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni2=1+value(m.ni2)
+ return m.dXNDdt[i] - (m.Q/m.Vw1)*(m.z42[i]-m.z2[i]) - m.rXND[i]==0
+ m.cz20=Constraint(m.t, rule=_cz20)
+
+
+
+ #organic nitrogen particulate biofilm
+ def _cz21(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni17=1+value(m.ni17)
+ return m.dXNDbdt[i] - m.rXNDb[i] ==0
+ m.cz21=Constraint(m.t, rule=_cz21)
+
+
+
+
+ # Inert soluble organic matter
+ def _cz22(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni28=1+value(m.ni28)
+ return m.dzSIdt2[i] - (m.Q/m.Vw1)*(m.z47[i]-m.z28[i]) ==0
+ m.cz22=Constraint(m.t, rule=_cz22)
+
+
+ # readily organic soluble
+ def _cz23(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni26=1+value(m.ni26)
+ return m.dzSdt2[i] - (m.Q/m.Vw1)*(m.z45[i]-m.z26[i]) - m.recs[i] ==0
+ m.cz23=Constraint(m.t, rule=_cz23)
+
+
+ #Inert particulate organic matter
+ def _cz24(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni29=1+value(m.ni29)
+ return m.dzXIdt2[i] - (m.Q/m.Vw1)*(0.52*((m.Q1+m.Q4[i])/m.Q1)*m.z48[i]-m.z29[i]) - m.recIX[i] ==0
+ m.cz24=Constraint(m.t, rule=_cz24)
+
+ #slowly organic matter
+ def _cz25(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni27=1+value(m.ni27)
+ return m.dzXdt2[i] - (m.Q/m.Vw1)*(0.52*((m.Q1+m.Q4[i])/m.Q1)*m.z46[i]-m.z27[i]) - m.recX[i] ==0
+ m.cz25=Constraint(m.t, rule=_cz25)
+
+
+ #Heterotrophic
+ def _cz26(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni30=1+value(m.ni30)
+ return m.dXHdt2[i] - (m.Q/m.Vw1)*(m.z52[i]-m.z30[i]) - m.recBH[i] ==0
+ m.cz26=Constraint(m.t, rule=_cz26)
+
+
+ #AOB
+ def _cz27(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni31=1+value(m.ni31)
+ return m.dAdt2[i] - (m.Q/m.Vw1)*(m.z51[i]-m.z31[i]) - m.recAOB[i] ==0
+ m.cz27=Constraint(m.t, rule=_cz27)
+
+
+ #NOB
+ def _cz28(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni65=1+value(m.ni65)
+ return m.ddAdt2[i] - (m.Q/m.Vw1)*(m.z66[i]-m.z65[i]) - m.recNOB[i] ==0
+ m.cz28=Constraint(m.t, rule=_cz28)
+
+
+ #Biomass decay
+ def _cz29(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni33=1+value(m.ni33)
+ return m.dXPdt2[i] - (m.Q/m.Vw1)*(0.52*((m.Q1+m.Q4[i])/m.Q1)*m.z50[i]-m.z33[i]) - m.recP[i] ==0
+ m.cz29=Constraint(m.t, rule=_cz29)
+
+
+ #Inert particulate organic matter biofilm
+ def _cz30(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni38=1+value(m.ni38)
+ return m.dzXIbdt2[i] - m.recIXb[i] ==0
+ m.cz30=Constraint(m.t, rule=_cz30)
+
+
+
+ #slowly organic matter biofilm
+ def _cz31(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni36=1+value(m.ni36)
+ return m.dzXbdt2[i] - m.recXb[i] ==0
+ m.cz31=Constraint(m.t, rule=_cz31)
+
+
+ #Heterotrophic biofilm
+ def _cz32(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni34=1+value(m.ni34)
+ return m.dXHbdt2[i] - m.recBHb[i] ==0
+ m.cz32=Constraint(m.t, rule=_cz32)
+
+
+ #AOB biofilm
+ def _cz33(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni35=1+value(m.ni35)
+ return m.dAbdt2[i] - m.recAOBb[i] ==0
+ m.cz33=Constraint(m.t, rule=_cz33)
+
+
+ #NOB biofilm
+ def _cz34(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni68=1+value(m.ni68)
+ return m.ddAbdt2[i] - m.recNOBb[i] ==0
+ m.cz34=Constraint(m.t, rule=_cz34)
+
+
+
+ #Biomass decay biofilm
+ def _cz35(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni39=1+value(m.ni39)
+ return m.dXPbdt2[i] - m.recPb[i] ==0
+ m.cz35=Constraint(m.t, rule=_cz35)
+
+
+ #NO3
+ def _cz37(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni32=1+value(m.ni32)
+ return m.dzNOdt2[i] - (m.Q/m.Vw1)*(m.z53[i] -m.z32[i]) - m.recNO3[i] ==0
+ m.cz37=Constraint(m.t, rule=_cz37)
+
+
+ #NO2
+ def _cz38(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni62=1+value(m.ni62)
+ return m.ddzNOdt2[i] - (m.Q/m.Vw1)*(m.z63[i] -m.z62[i]) - m.recNO2[i] ==0
+ m.cz38=Constraint(m.t, rule=_cz38)
+
+
+ #ammonia
+ def _cz39(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni20=1+value(m.ni20)
+ return m.dNHdt2[i] - (m.Q/m.Vw1)* (m.z40[i]-m.z20[i]) - m.recamm[i] ==0
+ m.cz39=Constraint(m.t, rule=_cz39)
+
+
+
+ #organic nitrogen soluble
+ def _cz40(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni21=1+value(m.ni21)
+ return m.dSNDdt2[i] - (m.Q/m.Vw1)* (m.z41[i] -m.z21[i]) - m.recSND[i] ==0
+ m.cz40=Constraint(m.t, rule=_cz40)
+
+
+
+ #organic nitrogen particulate
+ def _cz41(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni22=1+value(m.ni22)
+ return m.dXNDdt2[i] - (m.Q/m.Vw1)*(m.z42[i]-m.z22[i]) - m.recXND[i] ==0
+ m.cz41=Constraint(m.t, rule=_cz41)
+
+
+ #organic nitrogen particulate biofilm
+ def _cz42(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni37=1+value(m.ni37)
+ return m.dXNDbdt2[i] - m.recXNDb[i] ==0
+ m.cz42=Constraint(m.t, rule=_cz42)
+
+ def _weight(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni69=1+value(m.ni69)
+ return m.dWdt[i] - (1-m.agro)*m.bgro[i]*(m.tau*m.delta[i]*m.nu[i]*m.F[i]*1000/m.z70[i]) + m.Kgro * m.z69[i]**m.ngro == 0
+ m.weight=Constraint(m.t, rule=_weight)
+
+ def _coxy1(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni43=1+value(m.ni43)
+ return m.dy4dt[i] - ((m.Q1)/m.V)*(0.5*m.z23[i]+0.5*m.z3[i]) + ((m.Q1+m.Q4[i])/m.V)*(m.z43[i]) + (1/m.V)*(m.biom[i]*(3.05*10**(-6))* ((1.8* m.T + 32)**1.855)* (m.z69[i]/453.592)**(-0.138)) - m.kl[i] ==0 #O2
+ m.coxy1=Constraint(m.t, rule=_coxy1)
+
+
+ def _coxy2(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni3=1+value(m.ni3)
+ return m.dz02[i] - (m.Q/m.Vw1)*(m.z43[i]-m.z3[i]) - m.klBR1[i] - m.rO[i] ==0
+ m.coxy2=Constraint(m.t, rule=_coxy2)
+
+
+ def _coxy3(m,i):
+ if i==0:
+ return Constraint.Skip
+ else:
+ m.ni23=1+value(m.ni23)
+ return m.dz022[i] - (m.Q/m.Vw1)*(m.z43[i]-m.z23[i]) - m.klBR2[i] - m.recO[i] ==0
+ m.coxy3=Constraint(m.t, rule=_coxy3)
+
+
+ if flag==1:
+ disc=TransformationFactory('dae.finite_difference')
+ disc.apply_to(m,nfe=P,wrt=m.t,scheme='BACKWARD')
+ else:
+ disc=TransformationFactory('dae.collocation')
+ disc.apply_to(m,nfe=P*h,ncp=n_coll,wrt=m.t,scheme='LAGRANGE-RADAU')
+
+
+
+
+
+
+ def _init0(m):
+ return m.z0[0] == States[0]
+ m.init0 = Constraint(rule=_init0)
+
+ def _init1(m):
+ return m.z1[0] == States[1]
+ m.init1 = Constraint(rule=_init1)
+
+ def _init2(m):
+ return m.z2[0] == States[2]
+ m.init2 = Constraint(rule=_init2)
+
+ def _init6(m):
+ return m.z6[0] == States[6]
+ m.init6 = Constraint(rule=_init6)
+
+ def _init7(m):
+ return m.z7[0] == States[7]
+ m.init7 = Constraint(rule=_init7)
+
+ def _init8(m):
+ return m.z8[0] == States[8]
+ m.init8 = Constraint(rule=_init8)
+
+ def _init9(m):
+ return m.z9[0] == States[9]
+ m.init9 = Constraint(rule=_init9)
+
+ def _init10(m):
+ return m.z10[0] == States[10]
+ m.init10 = Constraint(rule=_init10)
+
+ def _init11(m):
+ return m.z11[0] == States[11]
+ m.init11 = Constraint(rule=_init11)
+
+ def _init12(m):
+ return m.z12[0] == States[12]
+ m.init12 = Constraint(rule=_init12)
+
+ def _init13(m):
+ return m.z13[0] == States[13]
+ m.init13 = Constraint(rule=_init13)
+
+ def _init14(m):
+ return m.z14[0] == States[14]
+ m.init14 = Constraint(rule=_init14)
+
+ def _init15(m):
+ return m.z15[0] == States[15]
+ m.init15 = Constraint(rule=_init15)
+
+ def _init16(m):
+ return m.z16[0] == States[16]
+ m.init16 = Constraint(rule=_init16)
+
+ def _init17(m):
+ return m.z17[0] == States[17]
+ m.init17 = Constraint(rule=_init17)
+
+ def _init18(m):
+ return m.z18[0] == States[18]
+ m.init18 = Constraint(rule=_init18)
+
+ def _init19(m):
+ return m.z19[0] == States[19]
+ m.init19 = Constraint(rule=_init19)
+
+ #Second reactor initial concentration
+
+ def _init20(m):
+ return m.z20[0] == States[20]
+ m.init20 = Constraint(rule=_init20)
+
+ def _init21(m):
+ return m.z21[0] == States[21]
+ m.init21 = Constraint(rule=_init21)
+
+ def _init22(m):
+ return m.z22[0] == States[22]
+ m.init22 = Constraint(rule=_init22)
+
+ def _init26(m):
+ return m.z26[0] == States[26]
+ m.init26 = Constraint(rule=_init26)
+
+ def _init27(m):
+ return m.z27[0] == States[27]
+ m.init27 = Constraint(rule=_init27)
+
+ def _init28(m):
+ return m.z28[0] == States[28]
+ m.init28 = Constraint(rule=_init28)
+
+ def _init29(m):
+ return m.z29[0] == States[29]
+ m.init29= Constraint(rule=_init29)
+
+ def _init30(m):
+ return m.z30[0] == States[30]
+ m.init30 = Constraint(rule=_init30)
+
+ def _init31(m):
+ return m.z31[0] == States[31]
+ m.init31 = Constraint(rule=_init31)
+
+ def _init32(m):
+ return m.z32[0] == States[32]
+ m.init32 = Constraint(rule=_init32)
+
+ def _init33(m):
+ return m.z33[0] == States[33]
+ m.init33 = Constraint(rule=_init33)
+
+ def _init34(m):
+ return m.z34[0] == States[34]
+ m.init34 = Constraint(rule=_init34)
+
+ def _init35(m):
+ return m.z35[0] == States[35]
+ m.init35 = Constraint(rule=_init35)
+
+ def _init36(m):
+ return m.z36[0] == States[36]
+ m.init36 = Constraint(rule=_init36)
+
+ def _init37(m):
+ return m.z37[0] == States[37]
+ m.init37 = Constraint(rule=_init37)
+
+ def _init38(m):
+ return m.z38[0] == States[38]
+ m.init38 = Constraint(rule=_init38)
+
+ def _init39(m):
+ return m.z39[0] == States[39]
+ m.init39 = Constraint(rule=_init39)
+
+ def _init40(m):
+ return m.z40[0] == States[40]
+ m.init40 = Constraint(rule=_init40)
+
+ def _init41(m):
+ return m.z41[0] == States[41]
+ m.init41 = Constraint(rule=_init41)
+
+ def _init42(m):
+ return m.z42[0] == States[42]
+ m.init42 = Constraint(rule=_init42)
+
+ def _init45(m):
+ return m.z45[0] == States[45]
+ m.init45 = Constraint(rule=_init45)
+
+ def _init46(m):
+ return m.z46[0] == States[46]
+ m.init46 = Constraint(rule=_init46)
+
+ def _init47(m):
+ return m.z47[0] == States[47]
+ m.init47 = Constraint(rule=_init47)
+
+ def _init48(m):
+ return m.z48[0] == States[48]
+ m.init48 = Constraint(rule=_init48)
+
+ def _init50(m):
+ return m.z50[0] == States[50]
+ m.init50 = Constraint(rule=_init50)
+
+ def _init51(m):
+ return m.z51[0] == States[51]
+ m.init51 = Constraint(rule=_init51)
+
+ def _init52(m):
+ return m.z52[0] == States[52]
+ m.init52 = Constraint(rule=_init52)
+
+ def _init53(m):
+ return m.z53[0] == States[53]
+ m.init53 = Constraint(rule=_init53)
+
+ def _init54(m):
+ return m.z54[0] == States[54]
+ m.init54 = Constraint(rule=_init54)
+
+ def _init55(m):
+ return m.z55[0] == States[55]
+ m.init55 = Constraint(rule=_init55)
+
+ def _init56(m):
+ return m.z56[0] == States[56]
+ m.init56 = Constraint(rule=_init56)
+
+ def _init57(m):
+ return m.z57[0] == States[57]
+ m.init57 = Constraint(rule=_init57)
+
+ def _init58(m):
+ return m.z58[0] == States[58]
+ m.init58 = Constraint(rule=_init58)
+
+ def _init59(m):
+ return m.z59[0] == States[59]
+ m.init59 = Constraint(rule=_init59)
+
+ def _init60(m):
+ return m.z60[0] == States[60]
+ m.init60 = Constraint(rule=_init60)
+
+ def _init61(m):
+ return m.z61[0] == States[61]
+ m.init61 = Constraint(rule=_init61)
+
+ def _init62(m):
+ return m.z62[0] == States[62]
+ m.init62 = Constraint(rule=_init62)
+
+ def _init63(m):
+ return m.z63[0] == States[63]
+ m.init63 = Constraint(rule=_init63)
+
+ def _init64(m):
+ return m.z64[0] == States[64]
+ m.init64 = Constraint(rule=_init64)
+
+ def _init65(m):
+ return m.z65[0] == States[65]
+ m.init65 = Constraint(rule=_init65)
+
+ def _init66(m):
+ return m.z66[0] == States[66]
+ m.init66 = Constraint(rule=_init66)
+
+ def _init67(m):
+ return m.z67[0] == States[67]
+ m.init67 = Constraint(rule=_init67)
+
+ def _init68(m):
+ return m.z68[0] == States[68]
+ m.init68 = Constraint(rule=_init68)
+
+ def _init69(m):
+ return m.z69[0] == States[69]
+ m.init69 = Constraint(rule=_init69)
+
+ def _init3(m):
+ return m.z3[0] == States[3]
+ m.init3 = Constraint(rule=_init3)
+
+ def _init23(m):
+ return m.z23[0] == States[23]
+ m.init23 = Constraint(rule=_init23)
+
+ def _init43(m):
+ return m.z43[0] == States[43]
+ m.init43 = Constraint(rule=_init43)
+
+
+ def _addcon_kl(m,i):
+ if i==0:
+ return m.kl[i]-inputs[1]==0
+ else:
+ return Constraint.Skip
+ m.addcon_kl=Constraint(m.t, rule=_addcon_kl)
+
+
+
+ def _addcon_klBR1(m,i):
+ if i==0:
+ return m.klBR1[i]-inputs[2]==0
+ else:
+ return Constraint.Skip
+ m.addcon_klBR1=Constraint(m.t, rule=_addcon_klBR1)
+
+
+
+ def _addcon_klBR2(m,i):
+ if i==0:
+ return m.klBR2[i]-inputs[4]==0
+ else:
+ return Constraint.Skip
+ m.addcon_klBR2=Constraint(m.t, rule=_addcon_klBR2)
+
+
+ def _addcon_Q4(m,i):
+ if i==0:
+ return m.Q4[i]-inputs[3]==0
+ else:
+ return Constraint.Skip
+ m.addcon_Q4=Constraint(m.t, rule=_addcon_Q4)
+
+
+
+
+
+ def _CVcons1(m,i):
+ return m.z3[i] >=0.007
+ m.CVcons1= Constraint(m.t,rule=_CVcons1)
+
+ def _CVcons2(m,i):
+ return m.z23[i] >=0.007
+ m.CVcons2= Constraint(m.t,rule=_CVcons2)
+
+ def _CVcons3(m,i):
+ return m.z43[i] >=0.007
+ m.CVcons3= Constraint(m.t,rule=_CVcons3)
+
+
+ def _CVcons4(m,i):
+ return m.z12[i] <=0.075
+ m.CVcons4= Constraint(m.t,rule=_CVcons4)
+
+ def _CVcons5(m,i):
+ return m.z32[i] <=0.075
+ m.CVcons5= Constraint(m.t,rule=_CVcons5)
+
+ def _CVcons6(m,i):
+ return m.z53[i] <=0.075
+ m.CVcons6= Constraint(m.t,rule=_CVcons6)
+
+
+
+ def _CVcons7(m,i):
+ return m.z0[i] <=0.013
+ m.CVcons7= Constraint(m.t,rule=_CVcons7)
+
+ def _CVcons8(m,i):
+ return m.z20[i] <=0.013
+ m.CVcons8= Constraint(m.t,rule=_CVcons8)
+
+ def _CVcons9(m,i):
+ return m.z40[i] <=0.013
+ m.CVcons9= Constraint(m.t,rule=_CVcons9)
+
+
+ m.obj=Objective(expr=(7.35-0.2)*m.z69[h*P]*m.z70[h*P]/1000 - sum(m.kl[i]*0.04*0.1 + 1.55*m.F[i]*0.1 for i in m.t) - sum( 0.0001*(m.kl[i]-m.kl[m.t.next(i)])**2 + 0.001*(m.F[i]-m.F[m.t.next(i)])**2 + 0.001*(m.Q4[i]-m.Q4[m.t.next(i)])**2 +0.0001*(m.klBR1[i]-m.klBR1[m.t.next(i)])**2+ 0.0001*(m.klBR2[i]-m.klBR2[m.t.next(i)])**2 for i in sorted(m.t)[:-1]) , sense = maximize)
+
+
+ return m
+
+
+
+
\ No newline at end of file
--
GitLab