"""
## GAMSSOURCE: https://www.gams.com/latest/noalib_ml/libhtml/noalib_protein.html
## LICENSETYPE: Requires license
## MODELTYPE: NLP
Optimal production of secreted protein in a fed-batch reactor.
Park, S., Ramirez, W.F., Optimal production of secreted protein in fed-batch
reactors. A.I.Ch.E. Journal, 34, 1988, pp.1550-1558.
"""
from __future__ import annotations
import os
import gamspy.math as gams_math
from gamspy import Alias
from gamspy import Container
from gamspy import Equation
from gamspy import Model
from gamspy import Options
from gamspy import Parameter
from gamspy import Set
from gamspy import Variable
def main():
m = Container(
system_directory=os.getenv("SYSTEM_DIRECTORY", None),
delayed_execution=int(os.getenv("DELAYED_EXECUTION", False)),
)
n = 500
# SET #
nh = Set(
m,
name="nh",
records=[str(i) for i in range(n + 1)],
description="Number of subintervals",
)
# ALIAS #
k = Alias(m, name="k", alias_with=nh)
# SCALARS #
tf = Parameter(m, name="tf", records=10, description="final time")
x1_0 = Parameter(
m, name="x1_0", records=1, description="initial value for x1"
)
x2_0 = Parameter(
m, name="x2_0", records=5, description="initial value for x2"
)
x3_0 = Parameter(
m, name="x3_0", records=0, description="initial value for x3"
)
x4_0 = Parameter(
m, name="x4_0", records=0, description="initial value for x4"
)
x5_0 = Parameter(
m, name="x5_0", records=1, description="initial value for x5"
)
h = Parameter(m, name="h")
h[...] = tf / n
# VARIABLES #
x1 = Variable(m, name="x1", domain=nh)
x2 = Variable(m, name="x2", domain=nh)
x3 = Variable(m, name="x3", domain=nh)
x4 = Variable(m, name="x4", domain=nh)
x5 = Variable(m, name="x5", domain=nh)
u = Variable(m, name="u", domain=nh, description="control variable")
a1 = Variable(m, name="a1", domain=nh)
a2 = Variable(m, name="a2", domain=nh)
a3 = Variable(m, name="a3", domain=nh)
# EQUATIONS #
state1 = Equation(
m,
name="state1",
type="regular",
domain=nh,
description="state equation 1",
)
state2 = Equation(
m,
name="state2",
type="regular",
domain=nh,
description="state equation 2",
)
state3 = Equation(
m,
name="state3",
type="regular",
domain=nh,
description="state equation 3",
)
state4 = Equation(
m,
name="state4",
type="regular",
domain=nh,
description="state equation 4",
)
state5 = Equation(
m,
name="state5",
type="regular",
domain=nh,
description="state equation 5",
)
ea1 = Equation(m, name="ea1", type="regular", domain=nh)
ea2 = Equation(m, name="ea2", type="regular", domain=nh)
ea3 = Equation(m, name="ea3", type="regular", domain=nh)
eobj = x4[str(n)] * x5[str(n)] # Objective function
state1[nh[k.lead(1)]] = x1[k.lead(1)] == (
x1[k]
+ (h / 2)
* (
a1[k] * x1[k]
- u[k] * x1[k] / x5[k]
+ a1[k.lead(1)] * x1[k.lead(1)]
- u[k.lead(1)] * x1[k.lead(1)] / x5[k.lead(1)]
)
)
state2[nh[k.lead(1)]] = x2[k.lead(1)] == (
x2[k]
+ (h / 2)
* (
-7.3 * a1[k] * x1[k]
- u[k] * (x2[k] - 20) / x5[k]
- 7.3 * a1[k.lead(1)] * x1[k.lead(1)]
- u[k.lead(1)] * (x2[k.lead(1)] - 20) / x5[k.lead(1)]
)
)
state3[nh[k.lead(1)]] = x3[k.lead(1)] == (
x3[k]
+ (h / 2)
* (
a2[k] * x1[k]
- u[k] * x3[k] / x5[k]
+ a2[k.lead(1)] * x1[k.lead(1)]
- u[k.lead(1)] * x3[k.lead(1)] / x5[k.lead(1)]
)
)
state4[nh[k.lead(1)]] = x4[k.lead(1)] == (
x4[k]
+ (h / 2)
* (
a3[k] * (x3[k] - x4[k])
- u[k] * x4[k] / x5[k]
+ a3[k.lead(1)] * (x3[k.lead(1)] - x4[k.lead(1)])
- u[k.lead(1)] * x4[k.lead(1)] / x5[k.lead(1)]
)
)
state5[nh[k.lead(1)]] = x5[k.lead(1)] == x5[k] + (h / 2) * (
u[k] + u[k.lead(1)]
)
ea1[nh[k]] = a1[k] == 21.87 * x2[k] / ((x2[k] + 0.4) * (x2[k] + 62.5))
ea2[nh[k]] = a2[k] == (x2[k] * gams_math.exp(-5 * x2[k])) / (0.1 + x2[k])
ea3[nh[k]] = a3[k] == 4.75 * a1[k] / (0.12 + a1[k])
# Initial point
x1.l[nh] = 1.0
x2.l[nh] = 5.0
x3.l[nh] = 0.0
x4.l[nh] = 0.0
x5.l[nh] = 1.0
u.l[nh] = 0.0
x1.fx["0"] = x1_0
x2.fx["0"] = x2_0
x3.fx["0"] = x3_0
x4.fx["0"] = x4_0
x5.fx["0"] = x5_0
# Bounds
u.lo[nh] = 0.0
u.up[nh] = 5
protein = Model(
m,
name="protein",
equations=m.getEquations(),
problem="nlp",
sense="max",
objective=eobj,
)
protein.solve(options=Options(time_limit=60000, iteration_limit=80000))
print(
"Objective Function Value: ", round(protein.objective_value, 4), "\n"
)
# REPORTING PARAMETER #
rep = Parameter(m, name="rep", domain=[nh, "*"])
rep[nh, "x1"] = x1.l[nh]
rep[nh, "x2"] = x2.l[nh]
rep[nh, "x3"] = x3.l[nh]
rep[nh, "x4"] = x4.l[nh]
rep[nh, "x5"] = x5.l[nh]
rep[nh, "u"] = u.l[nh]
rep.pivot().round(3).to_csv("protein_report.csv")
# End of protein
if __name__ == "__main__":
main()
