Skip to content
Snippets Groups Projects
Commit ff7f0ee6 authored by Gabriel Patron's avatar Gabriel Patron
Browse files

Upload New File

parent 84958ece
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MPC_RAS_tempdist_git.py 0 → 100644
+ 1995
− 0
View file @ ff7f0ee6
# 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment