Optimal control of coupled ODE model with economic cost function

[pavaninguva]

Hi! Amazing work with the code! It seems really capable but I seem to have some issues getting started and would be grateful for any help.

I have been trying optimize the inputs of the following bioreactor model,

$$ \begin{align} \frac{dT}{dt} &= \mu T - k_{1}V_{s}T + D(T_{in} - T) \nonumber \\ % \frac{dI_{s}}{dt} &= k_{1}V_{s}T -(k_{2}+D)I_{s} \nonumber \\ % \frac{dV_{s}}{dt} &= k_{3}I_{s} -k_{1}(T + I_{s})V_{s} -(k_{4}+D)V_{s} \nonumber \end{align} $$

The model has 3 states $T,I_{s},V_{s}$, 5 parameters $\mu, k_{1}-k_{4}$ and 2 decision variables / manipulated variables / inputs $D, T_{in}$. I would like to optimize the following economic cost function (corresponding to maximum bioreactor output.

$$ \max \int_{t=0}^{t=t_{end}} DV_{s} dt $$

I am basing my code off the examples in the documentation (Consumption Savings Problem and Pandemic Control)

using InfiniteOpt, Ipopt

#Define parameters and IC
x₀ = [5e6,0.0,1.25e5]
params = [2.12e-9,7.13e-3,168,0.035,0.027]
k₁,k₂,k₃,k₄,μ = params

#Set up optimizer and model
opt = Ipopt.Optimizer
ns = 500
m = InfiniteModel(opt)

#Specify variables 
@infinite_parameter(m,t in [0,200], num_supports=ns)

@variable(m,T>=0,Infinite(t))
@variable(m,Iₛ>=0, Infinite(t))
@variable(m,Vₛ>=0, Infinite(t))

@variable(m,0.0<=D<=0.1,Infinite(t))
@variable(m,0.0<=Tin<=1e7,Infinite(t))

#Objective function
@objective(m,Max,integral(D*Vₛ ,t))

#Initial conditions
@constraint(m,T(0) == x₀[1])
@constraint(m,Iₛ(0) == x₀[2])
@constraint(m,Vₛ(0) == x₀[3])

#ODE Equations
@constraint(m,c1,deriv(T,t) == μ*T -k₁*Vₛ*T +D*(Tin -T))
@constraint(m,c2,deriv(Iₛ,t) == k₁*Vₛ*T -(k₂+D)*Iₛ)
@constraint(m,c3, deriv(Vₛ,t) == k₃*Iₛ -k₁*(T+Iₛ)*Vₛ -(k₄+D)*Vₛ)

When I print the model i.e. print(m), I get a formulation that looks quite sensible:

Max ∫{t ∈ [0, 200]}[D(t)*Vₛ(t)]
Subject to
 T(t) ≥ 0.0, ∀ t ∈ [0, 200]
 Iₛ(t) ≥ 0.0, ∀ t ∈ [0, 200]
 Vₛ(t) ≥ 0.0, ∀ t ∈ [0, 200]
 D(t) ≥ 0.0, ∀ t ∈ [0, 200]
 D(t) ≤ 0.1, ∀ t ∈ [0, 200]
 Tin(t) ≥ 0.0, ∀ t ∈ [0, 200]
 Tin(t) ≤ 1.0e7, ∀ t ∈ [0, 200]
 T(0) ≥ 0.0
 T(0) = 5.0e6
 Iₛ(0) ≥ 0.0
 Iₛ(0) = 0.0
 Vₛ(0) ≥ 0.0
 Vₛ(0) = 125000.0
 c1 : 2.12e-9 T(t)*Vₛ(t) - Tin(t)*D(t) + D(t)*T(t) + ∂/∂t[T(t)] - 0.027 T(t) = 0.0, ∀ t ∈ [0, 200]
 c2 : -2.12e-9 T(t)*Vₛ(t) + D(t)*Iₛ(t) + ∂/∂t[Iₛ(t)] + 0.00713 Iₛ(t) = 0.0, ∀ t ∈ [0, 200]
 c3 : 2.12e-9 Vₛ(t)*T(t) + 2.12e-9 Vₛ(t)*Iₛ(t) + D(t)*Vₛ(t) + ∂/∂t[Vₛ(t)] - 168 Iₛ(t) + 0.035 Vₛ(t) = 0.0, ∀ t ∈ [0, 200]

However, when I run the optimizer, it seems to give some interesting behaviour. It would set $D$ to its maximum value and keep $T_{in}$ at zero which is very far from the optimal.

#Extract states
T_opt = value(T)
Is_opt = value(Iₛ)
Vs_opt = value(Vₛ)

#Extract inputs
D_opt = value(D)
Tin_opt = value(Tin)

ts = supports(t)
ix = 2:(length(ts)-1)

plot(ts[ix], D_opt[ix], lab="D")
plot!(ts[ix], Tin_opt[ix], lab="Tin")

I get the following output

Number of Iterations....: 438

                                   (scaled)                 (unscaled)
Objective...............:  -1.2257715258204984e+05    1.2257715258204984e+05
Dual infeasibility......:   4.7343146517503075e+04    4.7343146517503075e+04
Constraint violation....:   4.8752418489882615e+02    4.8752418489882615e+02
Variable bound violation:   9.9999966918229433e-09    9.9999966918229433e-09
Complementarity.........:   1.5026696464689725e-08    1.5026696464689725e-08
Overall NLP error.......:   5.3621089779661634e+02    4.7343146517503075e+04


Number of objective function evaluations             = 569
Number of objective gradient evaluations             = 109
Number of equality constraint evaluations            = 569
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 453
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 439
Total seconds in IPOPT                               = 9.900

EXIT: Converged to a point of local infeasibility. Problem may be infeasible.
QStandardPaths: XDG_RUNTIME_DIR not set, defaulting to '/tmp/runtime-pavan'

and the following plot for $D$ and $T_{in}$

Running the same optimization in other optimal control packages in Python e.g. Gekko gives me something much more sensible

So I have the following questions:

1. Is there a problem with my model formulation? It is surprising that the optimizer is forcing $T_{in}$ to zero even when its not being penalized with an input cost.
2. How does infiniteopt.jl know what are the manipulated variables and what are the states of the system
3. Are there ways to constrain the rate of change of the manipulated variables i.e. $| \frac{dT_{in}}{dt} | < M$ as often it is desirable to avoid sharp changes.

Thank you for the help!

[pulsipher]

Hi @pavaninguva and welcome to the forum! Thanks for trying out InfiniteOpt!

1. The short answer is your model might have a math error. Notice that the final Ipopt status indicates that "EXIT: Converged to a point of local infeasibility. Problem may be infeasible." This means that the optimizer was unable to find a feasible (i.e., valid solution). Thus, the values reported are essentially meaningless, since they do not satisfy the constraints. For future reference, we can use termination_status to check this (see https://infiniteopt.github.io/InfiniteOpt.jl/stable/guide/result/#Basic-Usage). Now the infeasible result either means the model is infeasible or the solver was unable to find a region of feasibility from the default initial guess (i.e., setting all variables to 0).
2. InfiniteOpt.jl is for general infinite-dimensional optimization problems is not tailored specifically for optimal control. Hence, like pyomo.dae there is no distinction between states and controls, the optimizer can figure out the degrees of freedom. More specific tools like Gekko do make such a distinction, but apart from some additional error checking and simulation capabilities, this information is ultimately not communicated to Ipopt when it solves the discretized problem.
3. Yep, you can simply pose that condition as a constraint directly using deriv(Tin, t) directly in the constraint. Another common trick in optimal control is to use a quadratic penalty term on the control variables in the objective to help smooth out the optimal policy. Note that when plotting control trajectories in InfiniteOpt, you will often see a kink originating at t=0. This is because the controls at t=0 can be any value since de don't affect the solution at all so the initial guess (defaulting to 0) is often used.

Now at this point you are probably wondering why Ipopt is giving you an infeasible result. I have a few ideas here:

• This model is not well-posed for floating point calculations since you have very large numbers (e.g., 1e7) multiplying into very small numbers (e.g., 1e-9). This can induce numerical instability. I recommend scaling the problem (e.g., by changing the units) to resolve this problem.
• Even with a scaled model, this might struggle to converge with the default initial guesses (provided by the start keyword argument). Providing a reasonable initial guess that is feasible can make a big difference. This is well illustrated in Minimum time manouevering problem + boundaries #219. You can find more information about specifying initial guesses at https://infiniteopt.github.io/InfiniteOpt.jl/stable/guide/result/#Basic-Usage.
• If all of the above doesn't work then the model might have a fundamental formulation error. I would try only having a few constraints at a time to isolate the culprit.

I hope that helps!