mintOC - User contributions [en]

Lotka Volterra fishing problem (Julia Neural Network solve)

2021-09-20T15:20:00Z

ClemensZeile:

This is an implementation of the [[Lotka Volterra fishing problem]] using Julia and a neural network by Julius Martensen, see https://github.com/AlCap23/IMPRS2021

The control function was relaxed, i.e. <math> u(t) \in [0,1] </math>.

<source lang="AMPL">
using LinearAlgebra
using Plots

using OrdinaryDiffEq
using Flux
using DiffEqFlux

# Example for optimal control

controller = Chain(Dense(1, 10, relu), Dense(10,10, relu), Dense(10,10, relu), Dense(10, 1, σ))

pnn, f = Flux.destructure(controller)

function lotka!(du, u, p, t)
c = f(p)([t])[1]
du[1] = (1-0.4*c)*u[1] - u[1]*u[2]
du[2] = u[1]*u[2] - (1+0.2*c)*u[2]
du[3] = (u[1] - 1)^2 + (u[2] -1)^2
end

u0 = Float32[0.5; 0.7; 0.0]
tspan = (0.0f0, 12.0f0)
Δt = 12.0f0/100

prob = ODEProblem(lotka!, u0, tspan, pnn)
solution = solve(prob, Tsit5(), saveat = Δt)
length(solution)

function loss(p)
s_ = solve(prob, Tsit5(), p = p)
l = sum(abs, s_[end,end])
return l, s_
end

loss(pnn)

ps = typeof(solution)[]
ls = eltype(p_init)[]

callback = function (p, l, pred)
display(l)
push!(ps, pred)
push!(ls, l)
l <= 1.35 && return true
return false
end

# Steer away from initial position and try to avoid local minima
res_1 = DiffEqFlux.sciml_train(loss, pnn, cb = callback, maxiters = 100)

anm = @animate for i in 1:1:length(ps)
pl_ = plot(ps[i], ylim = (0, 4), title = "Iteration $i : $(ls[i])")
plot!(ps[i].t, 1.34408f0*ones(length(ps[i].t)), color = :black, style = :dash, label = "Optimum")
pl_
end

gif(anm, joinpath(pwd(), "figures", "lotka_volterra_optimal_control.gif"), fps = 10)

# Plot the control value
s_ = solve(prob, Tsit5(), p = res_1.u, saveat = Δt)
us = f(res_1.u)(permutedims(s_.t))

p1 = plot(s_, legend = :bottomright, ylabel = "u(t)")
plot!(s_.t, 1.34408f0*ones(length(s_.t)), color = :black, style = :dash, label = "Optimum")

plot(
p1,
plot(s_.t, us', xlabel = "t", ylabel = "w(t)", label = nothing, xticks = 0:2:12),
layout = (2,1), link = :x
)

savefig(joinpath(pwd(), "figures", "lotka_volterra_optimal_control_input.png"))

plot(ls, title = "Objective", xlabel = "Iterations", label = nothing)
savefig(joinpath(pwd(), "figures", "lotka_volterra_optimal_control_objective.png"))

</source>

<gallery caption="Reference solution plot" widths="280px" heights="240px" perrow="1">
Image:Lotka volterra optimal control julia NN.gif| Animated solution process
</gallery>

[[Category: Julia/JuMP]]

File:Lotka volterra optimal control julia NN.gif

2021-09-20T15:17:37Z

ClemensZeile: Solution of the relaxed Lotka-Volterra-OCP with Julia and NN by https://github.com/AlCap23/IMPRS2021/blob/master/figures/lotka_volterra_optimal_control.gif

Solution of the relaxed Lotka-Volterra-OCP with Julia and NN by
https://github.com/AlCap23/IMPRS2021/blob/master/figures/lotka_volterra_optimal_control.gif

Lotka Volterra fishing problem

2021-09-20T15:16:03Z

ClemensZeile: /* Source Code */

{{Dimensions
|nd = 1
|nx = 3
|nw = 1
|nre = 3
}}The '''Lotka Volterra fishing problem''' looks for an optimal fishing strategy to be performed on a fixed time horizon to bring the biomasses of both predator as prey fish to a prescribed steady state. The problem was set up as a small-scale benchmark problem.
The well known [http://en.wikipedia.org/wiki/Lotka_volterra Lotka Volterra equations] for a predator-prey system have been augmented by an additional linear term, relating to fishing by man. The control can be regarded both in a relaxed, as in a discrete manner, corresponding to a part of the fleet, or the full fishing fleet.

The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. The interior point equality conditions fix the initial values of the differential states.

The optimal integer control functions shows [[:Category:Chattering|chattering]] behavior, making the Lotka Volterra fishing problem an ideal candidate for benchmarking of algorithms.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; c_0 x_0 \; w, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; c_1 x_1 \; w, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x(0) &=& (0.5, 0.7, 0)^T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. The decision, whether the fishing fleet is actually fishing at time <math>t</math> is denoted by <math>w(t)</math>.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_0, c_1) &=& (0.4, 0.2).
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of [http://en.wikipedia.org/wiki/Pontryagin%27s_minimum_principle Pontryagins maximum principle]. The optimal solution contains a singular arc, as can be seen in the plot of the optimal control. The two differential states and corresponding adjoint variables in the indirect approach are also displayed. A different approach to solving the relaxed problem is by using a direct method such as collocation or Bock's direct multiple shooting method. Optimal solutions for different control discretizations are also plotted in the leftmost figure.

The optimal objective value of this relaxed problem is <math>x_2(t_f) = 1.34408</math>. As follows from MIOC theory <bib id="Sager2011d" /> this is the best lower bound on the optimal value of the original problem with the integer restriction on the control function. In other words, this objective value can be approximated arbitrarily close, if the control only switches often enough between 0 and 1. As no optimal solution exists, two suboptimal ones are shown, one with only two switches and an objective function value of <math>x_2(t_f) = 1.38276</math>, and one with 56 switches and <math>x_2(t_f) = 1.34416</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:lotkaRelaxedControls.png| Optimal relaxed control determined by an indirect approach and by a direct approach on different control discretization grids.
Image:lotkaindirektStates.png| Differential states and corresponding adjoint variables in the indirect approach.
Image:lotka2Switches.png| Control and differential states with only two switches.
Image:lotka56Switches.png| Control and differential states with 56 switches.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:ACADO | ACADO code]] at [[Lotka Volterra fishing problem (ACADO)]]
* [[:Category:AMPL | AMPL code]] at [[Lotka Volterra fishing problem (AMPL)]]
* [[:Category:APMonitor | APMonitor code]] at [[Lotka Volterra fishing problem (APMonitor)]]
* [[:Category:Bocop | Bocop code]] at [[Lotka Volterra fishing problem (Bocop)]]
* [[:Category:Casadi | Casadi code]] at [[Lotka Volterra fishing problem (Casadi)]]
* [[:Category:Gekko | GEKKO Python code]] at [[Lotka Volterra fishing problem (GEKKO)]]
* [[:Category:JModelica | JModelica code]] at [[Lotka Volterra fishing problem (JModelica)]]
* [[:Category:Julia/JuMP | JuMP code]] at [[Lotka Volterra fishing problem (JuMP)]]
* [[:Category:Muscod | Muscod code]] at [[Lotka Volterra fishing problem (Muscod)]]
* [[:Category:switch | switch code]] at [[Lotka Volterra fishing problem (switch)]]
* [[:Category:TomDyn/PROPT | PROPT code]] at [[Lotka Volterra fishing problem (TomDyn/PROPT)]]
* [[:Category:Julia/JuMP | Julia code]] at [[Lotka Volterra fishing problem (Julia Neural Network solve)]]

== Variants ==

There are several alternative formulations and variants of the above problem, in particular

* a prescribed time grid for the control function <bib id="Sager2006" />, see also [[Lotka Volterra fishing problem (AMPL)]],
* a time-optimal formulation to get into a steady-state <bib id="Sager2005" />,
* the usage of a different target steady-state, as the one corresponding to <math> w(t) = 1</math> which is <math>(1 + c_1, 1 - c_0)</math>, see [[Lotka Volterra multi-arcs problem]]
* different fishing control functions for the two species, see [[Lotka Volterra Multimode fishing problem]]
* different fishing control functions that fish an absolute value from the two species, see [[Lotka Volterra absolute fishing problem]]
* a terminal constrained formulation, where a violation is penalized via slack variables, see [[Lotka Volterra (terminal constraint violation)]]
* different parameters and start values.

== Miscellaneous and Further Reading ==
The Lotka Volterra fishing problem was introduced by Sebastian Sager in a proceedings paper <bib id="Sager2006" /> and revisited in his PhD thesis <bib id="Sager2005" />. These are also the references to look for more details.

== References ==
<biblist />


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

Lotka Volterra fishing problem (Julia Neural Network solve)

2021-09-20T15:12:45Z

ClemensZeile: Created page with "This is an implementation of the Lotka Volterra fishing problem using Julia and a neural network by Julius Martensen, see https://github.com/AlCap23/IMPRS2021 The control..."

Three Tank multimode problem

2020-03-14T09:18:16Z

ClemensZeile:

{{Dimensions
|nd = 1
|nx = 3
|nw = 3
|nre = 2
}}This site describes a Double tank problem variant with three binary controls instead of only one control and three tanks, i.e., three differential states representing different compartments.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by

<math>
\begin{array}{llll}
\displaystyle \min_{x,w} & \displaystyle \int_{0}^{T} & k_1(x_2-k_2)^2 + k_3(x_3-k_4)^2 \; \text{d}t\\[1.5ex]
\mbox{s.t.} & \dot{x}_1 & = -\sqrt{x_1}+c_1 w_1 + c_2 w_2 - w_3 \sqrt{c_3 x_3}, \\[1.5ex]
& \dot{x}_2 & = \sqrt{x_1}-\sqrt{x_2}, \\[1.5ex]
& \dot{x}_3 & = \sqrt{x_2}-\sqrt{x_3}+w_3 \sqrt{c_3 x_3}, \\[1.5ex]
& x(0) & = (2,2,2)^T, \\[1.5ex]
& 1 & = \sum\limits_{i=1}^{3}w_i(t), \\
& w_i(t) &\in \{0, 1\}, \quad i=1\ldots 3.
\end{array}
</math>

== Parameters ==

These fixed values are used within the model.

<math> T=12, c_1=1, c_2=2, c_3=0.8, k_1=2, k_2=3, k_3 = 1, k_4 = 3. </math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of direct optimal control and the CIA decomposition. We denote the relaxed control values with <math>a(t)\in[0, 1]</math>.

The optimal objective value of the relaxed problem with <math> n_t=100, \, n_u=100 </math> is <math>8.775979</math>. The objective value of the binary controls obtained by Combinatorial Integral Approimation (CIA) is <math>8.789487</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:three_tank_relaxed_solution.png| Optimal relaxed controls and states determined by an direct approach with python 3.6 and CasADi, applied Multiple Shooting, 4th order Runge Kutta scheme and 100 discretization intervals.
Image:three_tank_binary_solution.png| According optimal binary controls and states determined by the direct approach. The relaxed controls were approximated by Combinatorial Integral Approximation.
Image:three_tank_rounding_solution.png| Binary and relaxed control values as part of the Combinatorial Integral Approximation problem
</gallery>

== Source Code ==

Model description is available in (using pycombina for solving the (CIA) rounding problem step):
* [[:Category:Casadi | Casadi code]] at [[Three Tank multimode problem (python/casadi)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]

Three Tank multimode problem (python/casadi)

2020-03-14T09:17:46Z

ClemensZeile: Created page with "=== Casadi === The model in Python code for a fixed control discretization grid using direct multiple shooting. We use pycombina for solving the (CIA) problem. <source lang="..."

=== Casadi ===

The model in Python code for a fixed control discretization grid using direct multiple shooting. We use pycombina for solving the (CIA) problem.
<source lang="Python">

# Three tank optimal control
# ----------------------
# An optimal control problem (OCP),
# solved with direct multiple-shooting.
#
# For more information on CasADi see: http://labs.casadi.org/OCP
from casadi import *
import pylab as pl
from pycombina import BinApprox, CombinaBnB
import matplotlib.pyplot as plt

############# Input variables #############
T = 12
N = 100 # number of control intervals

############# End of Input variables #############

opti = Opti() # Optimization problem

#------- parameter values --------------
# Model parameters
k_1 = 2
k_2 = 3
k_3 = 1
k_4 = 3

c_1 = 1
c_2 = 2
c_3 = 0.8
T = 12

# ---- decision variables ---------
X = opti.variable(4,N+1) # state trajectory
x1 = X[0,:]
x2 = X[1,:]
x3 = X[2,:]
x4 = X[3,:]

U = opti.variable(3,N) # control trajectory (throttle)
u1 = U[0,:]
u2 = U[1,:]
u3 = U[2,:]

# ---- objective ---------
opti.minimize(x4[N]) # Objective as Mayer term

# ---- dynamic constraints -------
f = lambda x,u: vertcat(c_1*u[0]+c_2*u[1]-(x[0])**0.5- u[2]* (c_3*x[0])**0.5 , (x[0])**0.5- (x[1])**0.5, (x[1])**0.5- (x[2])**0.5 +u[2]* (c_3*x[0])**0.5 , u[1]*0.0+ k_1*(x[1]-k_2)**2 + k_3*(x[2]-k_4)**2 ) # dx/dt = f(x,u)

dt = T/N # length of a control interval
for k in range(N): # loop over control intervals
# Runge-Kutta 4 integration
k1 = f(X[:,k], U[:,k])
k2 = f(X[:,k]+dt/2*k1, U[:,k])
k3 = f(X[:,k]+dt/2*k2, U[:,k])
k4 = f(X[:,k]+dt*k3, U[:,k])
x_next = X[:,k] + dt/6*(k1+2*k2+2*k3+k4)
opti.subject_to(X[:,k+1]==x_next) # close the gaps

# ---- path constraints -----------

opti.subject_to(opti.bounded(0,U,1)) # control is limited
opti.subject_to(u1+u2+u3 == 1) # SOS1 constraint

# ---- boundary conditions --------
opti.subject_to(x1[0]==2) # intitial values
opti.subject_to(x2[0]==2) # ...
opti.subject_to(x3[0]==2) # ...
opti.subject_to(x4[0]==0) # ... mayer term auxiliary state

#opti.subject_to(pos[-1]==1) # finish line at position 1

# ---- misc. constraints ----------
#opti.subject_to(T>=0) # Time must be positive

# ---- initial values for solver ---
opti.set_initial(x1, 2)
opti.set_initial(x2, 2)
opti.set_initial(x3, 2)
opti.set_initial(x4, 0)
opti.set_initial(u1, 1)
opti.set_initial(u2, 0)
opti.set_initial(u3, 0)

# ---- solve NLP ------
opti.solver("ipopt") # set numerical backend
sol = opti.solve() # actual solve

# ---- post-processing ------
from pylab import plot, step, figure, legend, show, spy

plt.plot(pl.linspace(0, 12, num=N+1),sol.value(x1),label="x1")
plt.plot(pl.linspace(0, 12, num=N+1),sol.value(x2),label="x2")
plt.plot(pl.linspace(0, 12, num=N+1),sol.value(x3),label="x3")

#plot(limit(sol.value(pos)),'r--',label="speed limit")
plt.step(pl.linspace(0, 12, num=N),sol.value(U[0,:]),'k',label="a1")
plt.step(pl.linspace(0, 12, num=N),sol.value(U[1,:]),label="a2")
plt.step(pl.linspace(0, 12, num=N),sol.value(U[2,:]),label="a3")

plt.xlabel('t')
plt.ylabel('state and control values')
plt.title('Differential state trajectories for relaxed controls')

#xlabel=("t")
plt.legend(loc="upper left")
#title="
plt.savefig('three_tank_relaxed_solution.png')

#################################################

### Setting up CIA problem

################################################

t = np.linspace(0,T,N+1)
b_rel = np.array([[min(1,max(0,x)) for x in sol.value(U[0,:])], [min(1,max(0,x)) for x in sol.value(U[1,:])], [min(1,max(0,x)) for x in sol.value(U[2,:])]])

binapprox = BinApprox(t = t, b_rel = b_rel, binary_threshold = 1e-3, \
off_state_included = False)

#binapprox.set_n_max_switches(n_max_switches = max_switches)
#binapprox.set_min_up_times(min_up_times = [min_up_time])
#binapprox.set_min_down_times(min_down_times = [min_down_time])

combina = CombinaBnB(binapprox)
combina.solve(use_warm_start=False, bnb_search_strategy='dfs', bnb_opts={'max_iter': 1000000})

b_bin_orig = pl.asarray(binapprox.b_bin)

#################################################

### Re-run Optimization Problem with fixed binary controls
### 'Dirty hack': just copy previous problem set up due to reuse issues of constraints

sol1_x1 = sol.value(x1)
sol1_x2 = sol.value(x2)
sol1_x3 = sol.value(x3)

opti2 = Opti() # Optimization problem

# ---- decision variables ---------
X = opti2.variable(4,N+1) # state trajectory
x1 = X[0,:]
x2 = X[1,:]
x3 = X[2,:]
x4 = X[3,:]

U = opti2.variable(3,N) # control trajectory (throttle)
u1 = U[0,:]
u2 = U[1,:]
u3 = U[2,:]

# ---- objective ---------
opti2.minimize(x4[N]) # Objective as Mayer term

# ---- dynamic constraints --------
f = lambda x,u: vertcat(c_1*u[0]+c_2*u[1]-(x[0])**0.5- u[2]* (c_3*x[0])**0.5 , (x[0])**0.5- (x[1])**0.5, (x[1])**0.5- (x[2])**0.5 +u[2]* (c_3*x[0])**0.5 , u[1]*0.0+ k_1*(x[1]-k_2)**2 + k_3*(x[2]-k_4)**2)

dt = T/N # length of a control interval
for k in range(N): # loop over control intervals
# Runge-Kutta 4 integration
k1 = f(X[:,k], U[:,k])
k2 = f(X[:,k]+dt/2*k1, U[:,k])
k3 = f(X[:,k]+dt/2*k2, U[:,k])
k4 = f(X[:,k]+dt*k3, U[:,k])
x_next = X[:,k] + dt/6*(k1+2*k2+2*k3+k4)
opti2.subject_to(X[:,k+1]==x_next) # close the gaps

# ---- path constraints -----------
opti2.subject_to(U==b_bin_orig) # control is fixed

# ---- boundary conditions --------
opti2.subject_to(x1[0]==2) # intitial values
opti2.subject_to(x2[0]==2) # ...
opti2.subject_to(x3[0]==2) # ...
opti2.subject_to(x4[0]==0) # ... mayer term auxiliary state

#opti.subject_to(pos[-1]==1) # finish line at position 1

# ---- misc. constraints ----------
#opti.subject_to(T>=0) # Time must be positive

# ---- initial values for solver ---
opti2.set_initial(sol.value_variables())
# ---- initial values for solver ---
opti2.set_initial(x1, sol1_x1)
opti2.set_initial(x2, sol1_x2)
opti2.set_initial(x3, sol1_x3)
opti2.set_initial(x4, 0)
#opti.set_initial(u1, 0)
#opti.set_initial(u2, 0)
#opti.set_initial(u3, 1)

# ---- solve NLP ------
opti2.solver("ipopt") # set numerical backend
sol2 = opti2.solve() # actual solve

######################
#### Plot the binary feasible trajectory solution
#####################

figure()

plt.plot(pl.linspace(0, 12, num=N+1),sol2.value(x1),label="x1")
plt.plot(pl.linspace(0, 12, num=N+1),sol2.value(x2),label="x2")
plt.plot(pl.linspace(0, 12, num=N+1),sol2.value(x3),label="x3")

#plot(limit(sol.value(pos)),'r--',label="speed limit")
plt.step(pl.linspace(0, 12, num=N),sol2.value(U[0,:]),'k',label="w1")
plt.step(pl.linspace(0, 12, num=N),sol2.value(U[1,:]),label="w2")
plt.step(pl.linspace(0, 12, num=N),sol2.value(U[2,:]),label="w3")

plt.xlabel('t')
plt.ylabel('state and control values')
plt.title('Differential state trajectories for binary controls')

#xlabel=("t")
plt.legend(loc="upper left")
#title="
show()
plt.savefig('three_tank_binary_solution.png')

######################
#### Plot also the (CIA) rounding trajectories
#####################

f, (ax1, ax2, ax3) = pl.subplots(3, sharex = True)

ax1.step(t[:-1], b_bin_orig[0,:], label = "w", color = "C0", where = "post")
# ax1.step(t[:-1], b_bin_red[0,:], label = "b_bin_red", color = "C0", linestyle = "dashed", where = "post")
ax1.scatter(t[:-1], b_rel[0,:], label = "a", color = "C0", marker = "x")
ax1.legend(loc = "upper left")
ax1.set_ylabel("w_1")

ax2.step(t[:-1], b_bin_orig[1,:], label = "w", color = "C1", where = "post")
# ax2.step(t[:-1], b_bin_red[1,:], label = "b_bin_red", color = "C1", linestyle = "dashed", where = "post")
ax2.scatter(t[:-1], b_rel[1,:], label = "a", color = "C1", marker = "x")
ax2.legend(loc = "upper left")
ax2.set_ylabel("w_2")

ax3.step(t[:-1], b_bin_orig[2,:], label = "w", color = "C2", where = "post")
# ax3.step(t[:-1], b_bin_red[2,:], label = "b_bin_red", color = "C2", linestyle = "dashed", where = "post")
ax3.scatter(t[:-1], b_rel[2,:], label = "a", color = "C2", marker = "x")
ax3.legend(loc = "upper left")
ax3.set_ylabel("w_3")
ax3.set_xlabel("t")

show()
plt.savefig('three_tank_rounding_solution.png')
################################################

<source>

Three Tank multimode problem

2020-03-14T09:13:31Z

ClemensZeile: /* Reference Solutions */

{{Dimensions
|nd = 1
|nx = 3
|nw = 3
|nre = 2
}}This site describes a Double tank problem variant with three binary controls instead of only one control and three tanks, i.e., three differential states representing different compartments.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by

<math>
\begin{array}{llll}
\displaystyle \min_{x,w} & \displaystyle \int_{0}^{T} & k_1(x_2-k_2)^2 + k_3(x_3-k_4)^2 \; \text{d}t\\[1.5ex]
\mbox{s.t.} & \dot{x}_1 & = -\sqrt{x_1}+c_1 w_1 + c_2 w_2 - w_3 \sqrt{c_3 x_3}, \\[1.5ex]
& \dot{x}_2 & = \sqrt{x_1}-\sqrt{x_2}, \\[1.5ex]
& \dot{x}_3 & = \sqrt{x_2}-\sqrt{x_3}+w_3 \sqrt{c_3 x_3}, \\[1.5ex]
& x(0) & = (2,2,2)^T, \\[1.5ex]
& 1 & = \sum\limits_{i=1}^{3}w_i(t), \\
& w_i(t) &\in \{0, 1\}, \quad i=1\ldots 3.
\end{array}
</math>

== Parameters ==

These fixed values are used within the model.

<math> T=12, c_1=1, c_2=2, c_3=0.8, k_1=2, k_2=3, k_3 = 1, k_4 = 3. </math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of direct optimal control and the CIA decomposition. We denote the relaxed control values with <math>a(t)\in[0, 1]</math>.

The optimal objective value of the relaxed problem with <math> n_t=100, \, n_u=100 </math> is <math>8.775979</math>. The objective value of the binary controls obtained by Combinatorial Integral Approimation (CIA) is <math>8.789487</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:three_tank_relaxed_solution.png| Optimal relaxed controls and states determined by an direct approach with python 3.6 and CasADi, applied Multiple Shooting, 4th order Runge Kutta scheme and 100 discretization intervals.
Image:three_tank_binary_solution.png| According optimal binary controls and states determined by the direct approach. The relaxed controls were approximated by Combinatorial Integral Approximation.
Image:three_tank_rounding_solution.png| Binary and relaxed control values as part of the Combinatorial Integral Approximation problem
</gallery>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]

Three Tank multimode problem

2020-03-14T09:12:39Z

ClemensZeile: /* Reference Solutions */

{{Dimensions
|nd = 1
|nx = 3
|nw = 3
|nre = 2
}}This site describes a Double tank problem variant with three binary controls instead of only one control and three tanks, i.e., three differential states representing different compartments.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by

<math>
\begin{array}{llll}
\displaystyle \min_{x,w} & \displaystyle \int_{0}^{T} & k_1(x_2-k_2)^2 + k_3(x_3-k_4)^2 \; \text{d}t\\[1.5ex]
\mbox{s.t.} & \dot{x}_1 & = -\sqrt{x_1}+c_1 w_1 + c_2 w_2 - w_3 \sqrt{c_3 x_3}, \\[1.5ex]
& \dot{x}_2 & = \sqrt{x_1}-\sqrt{x_2}, \\[1.5ex]
& \dot{x}_3 & = \sqrt{x_2}-\sqrt{x_3}+w_3 \sqrt{c_3 x_3}, \\[1.5ex]
& x(0) & = (2,2,2)^T, \\[1.5ex]
& 1 & = \sum\limits_{i=1}^{3}w_i(t), \\
& w_i(t) &\in \{0, 1\}, \quad i=1\ldots 3.
\end{array}
</math>

== Parameters ==

These fixed values are used within the model.

<math> T=12, c_1=1, c_2=2, c_3=0.8, k_1=2, k_2=3, k_3 = 1, k_4 = 3. </math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of direct optimal control and the CIA decomposition. We denote the relaxed control values with <math>a(t)\in[0, 1]</math>.

The optimal objective value of the relaxed problem with <math> n_t=100, \, n_u=100 </math> is <math>8.775979</math>. The objective value of the binary controls obtained by Combinatorial Integral Approimation (CIA) is <math>8.789487</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:three_tank_relaxed_solution.png| Optimal relaxed controls and states determined by an direct approach with python 3.6 and CasADi, applied Multiple Shooting, 4th order Runge Kutta scheme and <math>n_u=100</math>.
Image:three_tank_binary_solution.png| According optimal binary controls and states determined by the direct approach. The relaxed controls were approximated by Combinatorial Integral Approximation.
Image:three_tank_rounding_solution.png| Binary and relaxed control values as part of the Combinatorial Integral Approximation problem
</gallery>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]

File:Three tank rounding solution.png

2020-03-14T09:05:26Z

ClemensZeile: Solution obtained by pycombina for the three tank OC problem

Solution obtained by pycombina for the three tank OC problem

File:Three tank relaxed solution.png

2020-03-14T09:04:40Z

ClemensZeile: Solution of three tank OC problem with relaxed controls

Solution of three tank OC problem with relaxed controls

File:Three tank binary solution.png

2020-03-14T09:00:40Z

ClemensZeile: Solution of three tank OC problem with binary controls fixed

Solution of three tank OC problem with binary controls fixed

File:Three tank binary solution.pdf

2020-03-14T08:58:21Z

ClemensZeile: Solution of three tank OC problem with binary controls fixed

Solution of three tank OC problem with binary controls fixed

Three Tank multimode problem

2020-03-14T07:48:51Z

ClemensZeile: Created page with "{{Dimensions |nd = 1 |nx = 3 |nw = 3 |nre = 2 }}This site describ..."

{{Dimensions
|nd = 1
|nx = 3
|nw = 3
|nre = 2
}}This site describes a Double tank problem variant with three binary controls instead of only one control and three tanks, i.e., three differential states representing different compartments.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by

<math>
\begin{array}{llll}
\displaystyle \min_{x,w} & \displaystyle \int_{0}^{T} & k_1(x_2-k_2)^2 + k_3(x_3-k_4)^2 \; \text{d}t\\[1.5ex]
\mbox{s.t.} & \dot{x}_1 & = -\sqrt{x_1}+c_1 w_1 + c_2 w_2 - w_3 \sqrt{c_3 x_3}, \\[1.5ex]
& \dot{x}_2 & = \sqrt{x_1}-\sqrt{x_2}, \\[1.5ex]
& \dot{x}_3 & = \sqrt{x_2}-\sqrt{x_3}+w_3 \sqrt{c_3 x_3}, \\[1.5ex]
& x(0) & = (2,2,2)^T, \\[1.5ex]
& 1 & = \sum\limits_{i=1}^{3}w_i(t), \\
& w_i(t) &\in \{0, 1\}, \quad i=1\ldots 3.
\end{array}
</math>

== Parameters ==

These fixed values are used within the model.

<math> T=12, c_1=1, c_2=2, c_3=0.8, k_1=2, k_2=3, k_3 = 1, k_4 = 3. </math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of direct optimal control.

The optimal objective value of the relaxed problem with <math> n_t=12000, \, n_u=100 </math> is <math>2.59106823</math>. The objective value of the binary controls obtained by Combinatorial Integral Approimation (CIA) is <math>2.59121008</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:MmdoubletankRelaxed 12000 120 1.png| Optimal relaxed controls and states determined by an direct approach with ampl_mintoc (Radau collocation) and <math>n_t=12000, \, n_u=100</math>.
Image:MmdoubletankCIA_12000_120_1.png| Optimal binary controls and states determined by an direct approach (Radau collocation) with ampl_mintoc and <math>n_t=12000, \, n_u=100</math>. The relaxed controls were approximated by Combinatorial Integral Approximation.
</gallery>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]

Double Tank multimode problem

2020-03-14T07:37:53Z

ClemensZeile:

{{Dimensions
|nd = 1
|nx = 2
|nw = 3
|nre = 2
}}This site describes a Double tank problem variant with three binary controls instead of only one control.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by

<math>
\begin{array}{llll}
\displaystyle \min_{x,w} & \displaystyle \int_{0}^{T} & k_1(x_2-k_2)^2 \; \text{d}t\\[1.5ex]
\mbox{s.t.} & \dot{x}_1 & = \sum\limits_{i=1}^{3} c_{i}\; w_i,-\sqrt{x_1}, \\[1.5ex]
& \dot{x}_2 & = \sqrt{x_1}-\sqrt{x_2}, \\[1.5ex]
& x(0) & = (2,2)^T, \\[1.5ex]
& 1 & = \sum\limits_{i=1}^{3}w_i(t), \\
& w_i(t) &\in \{0, 1\}, \quad i=1\ldots 3.
\end{array}
</math>

== Parameters ==

These fixed values are used within the model.

<math> T=10, c_1=1, c_2=0.5, c_3=2, k_1=2, k_2=3. </math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of direct optimal control.

The optimal objective value of the relaxed problem with <math> n_t=12000, \, n_u=100 </math> is <math>2.59106823</math>. The objective value of the binary controls obtained by Combinatorial Integral Approimation (CIA) is <math>2.59121008</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:MmdoubletankRelaxed 12000 120 1.png| Optimal relaxed controls and states determined by an direct approach with ampl_mintoc (Radau collocation) and <math>n_t=12000, \, n_u=100</math>.
Image:MmdoubletankCIA_12000_120_1.png| Optimal binary controls and states determined by an direct approach (Radau collocation) with ampl_mintoc and <math>n_t=12000, \, n_u=100</math>. The relaxed controls were approximated by Combinatorial Integral Approximation.
</gallery>

== Source Code ==

Model description is available in
* [[:Category:AMPL | AMPL code]] at [[Double Tank multimode problem (AMPL)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]

Double Tank

2020-03-14T07:34:17Z

ClemensZeile: /* Variants */

{{Dimensions
|nd = 1
|nx = 2
|nw = 1
|nre = 2
}}

The double tank problem is a basic example for a switching system. It contains the dynamics of an upper and a lower tank, connected to each other with a pipe. The goal is to minimize the deviation of a certain fluid level <math>k_2</math> in the lower tank. The problem was introduced and discussed in a variety of publications for the optimal control of constrained switched systems, e.g. <bib id="Henrion2014" /> (Link: [http://homepages.laas.fr/henrion/papers/switch.pdf Henrion et al.]) and <bib id="Vasudevan2013" /> (Link: [http://epubs.siam.org/doi/pdf/10.1137/120901507 Vasudevan et al.]).

== Mathematical formulation ==

<math>
\begin{array}{llll}
\displaystyle \min_{\sigma} & \displaystyle \int_{0}^{T} & k_1(x_2-k_2)^2 \; \text{d}t\\[1.5ex]
\mbox{s.t.} & \dot{x}_1 & = c_{\sigma(t)}-\sqrt{x_1}, \\[1.5ex]
& \dot{x}_2 & = \sqrt{x_1}-\sqrt{x_2}, \\[1.5ex]
& x(0) & = (2,2)^T, \\[1.5ex]
& \sigma(t) & \in \{1,2\} \qquad &\text{for } t\in[0,T].\\[1.5ex]
\end{array}
</math>

The two states of the system correspond to the fluid levels of an upper and a lower tank.
The output of the upper tank flows into the lower tank, the output of the lower
tank exits the system, and the flow into the upper tank is restricted to either <math>c_1</math> [lt/s]
or <math>c_2</math> [lt/s]. The dynamics in each mode are then derived using [https://en.wikipedia.org/wiki/Torricelli's_law Torricelli’s law], as
shown in constraints 1 and 2. The objective of the optimization is to have the fluid level in the
lower tank equal to <math>k_2</math> [m], as reflected in the cost function.

== Parameters ==
In an exemplary test, the parameters were chosen to be:

{| class="wikitable"
|+Parameters
|-
|Symbol
|Value
|-
|<math>k_1</math>
|<math>2</math>
|-
|<math>k_2</math>
|<math>3</math>
|-
|<math>c_1</math>
|<math>1</math>
|-
|<math>c_2</math>
|<math>2</math>
|-
|<math>T</math>
|<math>10</math>
|}

== Reference Solution ==
The problem is solved in Python with [[:Category: Gekko | GEKKO]] or Matlab (version 2014b) with [[:Category: switch | Switch]]. The switch codes rely on [http://homepages.laas.fr/henrion/software/gloptipoly3/ GloptiPoly 3] and the solver [http://sedumi.ie.lehigh.edu/ SeDuMi 1.3] (optimization over symmetric cones).
By introducing the lifts <math>l_i=\sqrt{x_i}</math>, algebraically constrained as <math>l_i^2=x_i, \; l_i\geq 0,</math> the problem is recast with polynomial data. In this way way switch in connection with GloptiPoly3 can be applied. GloptiPoly 3 is a Matlab package developed by Didier Henrion for generalized problems of moments and needs polynomial input data. The exact code used to solve the problem can be found under [[Double Tank (switch)]]. The calculated objective is 4.7296 with the following trajectories of states and controls:
<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:State_sequence_Tank.jpg| State trajectory over time.
Image:Control_sequence_Tank.jpg| Control sequence over time.
</gallery>

== Variants ==

There are several alternative formulations and variants of the above problem, in particular

* different control functions for the two tanks, see [[Double Tank multimode problem]]
* a three tank problem with 3 modes, see [[Three Tank multimode problem]]

== Source Code ==

Model descriptions are available in

* [[:Category:Gekko | GEKKO Python code]] at [[Double Tank (GEKKO)]]
* [[:Category: switch | switch code]] at [[Double Tank (switch)]]

== References ==
<biblist />

CLAEYS, Mathieu; DAAFOUZ Jamal; HENRION Didier Modal occupation measures and LMI relaxations for nonlinear switched systems control. arXiv preprint arXiv:1404.4699 (2014)

VASUDEVAN, Ramanarayan, et al. Consistent Approximations for the Optimal Control of Constrained Switched Systems---Part 2: An Implementable Algorithm. SIAM Journal on Control and Optimization, 2013, 51. Jg., Nr. 6, S. 4484-4503.


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category: Tracking objective]]
[[Category: State dependent switches]]
[[Category:Bang bang]]
[[Category:Path-constrained arcs]]

Quadrotor (binary variant)

2019-10-14T14:29:39Z

ClemensZeile: Created page with "{{Dimensions |nd = 1 |nx = 6 |nw = 4 |nre = 6 }}This site describ..."

{{Dimensions
|nd = 1
|nx = 6
|nw = 4
|nre = 6
}}This site describes a Quadrotor helicoptor problem variant where the continuous control is replaced via outer convexification with binary controls.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x,u, w} & 5(x_1(t_f)-6)^2&+&5(x_3(t_f)-1)^2+(\sin(x_5(t_f)0.5))^2 +\int\limits_{t_0}^{t_f} 5( (w_2(\tau)+w_4(\tau)+w_6(\tau))^2 \ d \tau \\[1.5ex]
\mbox{s.t.}
& \dot{x}_1 & = & x_2(t), \\
& \dot{x}_2 & = & g \sin( x_5(t)) + \sum\limits_{i\in [4]}c_{1,i}w_i(t)\frac{\sin(x_5(t))}{M}, \\
& \dot{x}_3 & = & x_4(t), \\
& \dot{x}_4 & = & g \cos( x_5(t))-g+ \sum\limits_{i\in [4]}c_{1,i}w_i(t)\frac{\cos(x_5(t))}{M}, \\
& \dot{x}_5 & = & x_6(t), \\
& \dot{x}_6 & = & \sum\limits_{i\in [4]}c_{2,i}w_i(t)L \frac{1}{I} \\[1.5ex]
& x(0) &=& (0, 0, 1, 0 , 0, 0)^T, \\
& w_i(t) &\in& \{0, 1\}, i=1,\ldots,4 \\
& \sum\limits_{i=1}^{4}w_i(t) &=& 1, \\
& x_3(t) & \geq & 0, \quad t\in[t_0,t_f].
\end{array}
</math>


== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 7.5],\\
(g, M, L, I) &=& (9.8, 1.3, 0.305, 0.0605),\\
c_1 &=& (0,0.001,0,0),\\
c_2 &=& (0,0,-0.001,0.001),
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> is in the continuous interval <math>[0, 1]</math> rather than being binary, the optimal solution can be determined by means of direct optimal control.

The optimal objective value of the relaxed problem with <math> n_t=12000, \, n_u=25 </math> is <math>13.0907346</math>. The objective value of the solution with binary controls obtained by Combinatorial Integral Approximation (CIA) is <math>15.5787932</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:Quadrotor_bin_Relaxed_12000_25.pdf| Optimal relaxed controls and states determined by an direct approach with ampl_mintoc (Radau collocation) and <math>n_t=12000, \, n_u=25</math>.
Image:Quadrotor_bin_CIA_12000_25.pdf| Differential states and binary cotnrol determined by an direct approach (Radau collocation) with ampl_mintoc and <math>n_t=12000, \, n_u=25</math>. The relaxed controls were approximated by Combinatorial Integral Approximation.
</gallery>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

File:Quadrotor bin CIA 12000 25.pdf

2019-10-14T14:28:48Z

ClemensZeile: Numerical results for the outer convexified problem "Quadrotor Helicoptor" with nt=12000, nu=25 and binary controls obtained via CIA decomposition.

Numerical results for the outer convexified problem "Quadrotor Helicoptor" with nt=12000, nu=25 and binary controls obtained via CIA decomposition.

File:Quadrotor bin Relaxed 12000 25.pdf

2019-10-14T14:28:16Z

ClemensZeile: Numerical results for the outer convexified problem "Quadrotor Helicoptor" with nt=12000, nu=25 and relaxed binary controls.

Numerical results for the outer convexified problem "Quadrotor Helicoptor" with nt=12000, nu=25 and relaxed binary controls.

Quadrotor helicopter control problem

2019-10-14T13:44:09Z

ClemensZeile:

{{Dimensions
|nd = 1
|nx = 6
|nu = 1
|nw = 3
|nre = 6
|nrc = 3
}}The mixed-integer optimal control problem of a quadrotor helicopter in two dimensions is taken from (Link: [https://pdfs.semanticscholar.org/75a0/211476ddc21363cfb3262c04d18794ad06ef.pdf Gillula et al.]) and from (Link: [http://epubs.siam.org/doi/pdf/10.1137/120901507 Vasudevan et al.]). The evolution
of the quadrotor can be defined with respect to a fixed two dimensional reference
frame using six dimensions, where the first three dimensions represent the position
along a horizontal axis, the position along the vertical axis, and the roll angle of the
helicopter, respectively, and the last three dimensions represent the time derivative of
the first three dimensions.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x,u, w} & 5(x_1(t_f)-6)^2&+&5(x_3(t_f)-1)^2+(\sin(x_5(t_f)0.5))^2 +\int\limits_{t_0}^{t_f} 5u(\tau)^2 \ d \tau \\[1.5ex]
\mbox{s.t.}
& \dot{x}_1 & = & x_2(t), \\
& \dot{x}_2 & = & g \sin( x_5(t)) + w_1(t)u(t)\frac{\sin(x_5(t))}{M}, \\
& \dot{x}_3 & = & x_4(t), \\
& \dot{x}_4 & = & g \cos( x_5(t))-g+ w_1(t)u(t)\frac{\cos(x_5(t))}{M}, \\
& \dot{x}_5 & = & x_6(t), \\
& \dot{x}_6 & = & -w_2(t)L \frac{u(t)}{I}+w_3(t)L \frac{u(t)}{I} \\[1.5ex]
& x(0) &=& (0, 0, 1, 0 , 0, 0)^T, \\
& w_i(t) &\in& \{0, 1\}, i=1,\ldots,3 \\
& \sum\limits_{i=1}^{3}w_i(t) &=& 1, \\
& u(t) & \in& [0,0.001], \quad t\in[t_0,t_f],\\
& x_3(t) & \geq & 0, \quad t\in[t_0,t_f].
\end{array}
</math>


== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 7.5],\\
(g, M, L, I) &=& (9.8, 1.3, 0.305, 0.0605),
\end{array}
</math>

== Reference Solutions ==

A reference solution can be found in [http://epubs.siam.org/doi/pdf/10.1137/120901507 Vasudevan et al.] based on the embedding transformation technique for switched systems.

== Variants ==

There are several alternative formulations and variants of the above problem, in particular

* [[Quadrotor (binary variant)]]: The quadrotor helicoptor problem, where the continuous control is replaced via partial outer convexification by binary controls.


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]

Quadrotor helicopter control problem

2019-10-14T13:43:36Z

ClemensZeile:

{{Dimensions
|nd = 1
|nx = 6
|nu = 1
|nw = 3
|nre = 6
|nrc = 3
}}The mixed-integer optimal control problem of a quadrotor helicopter in two dimensions is taken from (Link: [https://pdfs.semanticscholar.org/75a0/211476ddc21363cfb3262c04d18794ad06ef.pdf Gillula et al.]).
and from (Link: [http://epubs.siam.org/doi/pdf/10.1137/120901507 Vasudevan et al.]). The evolution
of the quadrotor can be defined with respect to a fixed two dimensional reference
frame using six dimensions, where the first three dimensions represent the position
along a horizontal axis, the position along the vertical axis, and the roll angle of the
helicopter, respectively, and the last three dimensions represent the time derivative of
the first three dimensions.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x,u, w} & 5(x_1(t_f)-6)^2&+&5(x_3(t_f)-1)^2+(\sin(x_5(t_f)0.5))^2 +\int\limits_{t_0}^{t_f} 5u(\tau)^2 \ d \tau \\[1.5ex]
\mbox{s.t.}
& \dot{x}_1 & = & x_2(t), \\
& \dot{x}_2 & = & g \sin( x_5(t)) + w_1(t)u(t)\frac{\sin(x_5(t))}{M}, \\
& \dot{x}_3 & = & x_4(t), \\
& \dot{x}_4 & = & g \cos( x_5(t))-g+ w_1(t)u(t)\frac{\cos(x_5(t))}{M}, \\
& \dot{x}_5 & = & x_6(t), \\
& \dot{x}_6 & = & -w_2(t)L \frac{u(t)}{I}+w_3(t)L \frac{u(t)}{I} \\[1.5ex]
& x(0) &=& (0, 0, 1, 0 , 0, 0)^T, \\
& w_i(t) &\in& \{0, 1\}, i=1,\ldots,3 \\
& \sum\limits_{i=1}^{3}w_i(t) &=& 1, \\
& u(t) & \in& [0,0.001], \quad t\in[t_0,t_f],\\
& x_3(t) & \geq & 0, \quad t\in[t_0,t_f].
\end{array}
</math>


== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 7.5],\\
(g, M, L, I) &=& (9.8, 1.3, 0.305, 0.0605),
\end{array}
</math>

== Reference Solutions ==

A reference solution can be found in [http://epubs.siam.org/doi/pdf/10.1137/120901507 Vasudevan et al.] based on the embedding transformation technique for switched systems.

== Variants ==

There are several alternative formulations and variants of the above problem, in particular

* [[Quadrotor (binary variant)]]: The quadrotor helicoptor problem, where the continuous control is replaced via partial outer convexification by binary controls.


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]

Lotka Volterra fishing problem

2019-10-14T13:42:03Z

ClemensZeile: /* Variants */

{{Dimensions
|nd = 1
|nx = 3
|nw = 1
|nre = 3
}}The '''Lotka Volterra fishing problem''' looks for an optimal fishing strategy to be performed on a fixed time horizon to bring the biomasses of both predator as prey fish to a prescribed steady state. The problem was set up as a small-scale benchmark problem.
The well known [http://en.wikipedia.org/wiki/Lotka_volterra Lotka Volterra equations] for a predator-prey system have been augmented by an additional linear term, relating to fishing by man. The control can be regarded both in a relaxed, as in a discrete manner, corresponding to a part of the fleet, or the full fishing fleet.

The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. The interior point equality conditions fix the initial values of the differential states.

The optimal integer control functions shows [[:Category:Chattering|chattering]] behavior, making the Lotka Volterra fishing problem an ideal candidate for benchmarking of algorithms.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; c_0 x_0 \; w, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; c_1 x_1 \; w, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x(0) &=& (0.5, 0.7, 0)^T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. The decision, whether the fishing fleet is actually fishing at time <math>t</math> is denoted by <math>w(t)</math>.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_0, c_1) &=& (0.4, 0.2).
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of [http://en.wikipedia.org/wiki/Pontryagin%27s_minimum_principle Pontryagins maximum principle]. The optimal solution contains a singular arc, as can be seen in the plot of the optimal control. The two differential states and corresponding adjoint variables in the indirect approach are also displayed. A different approach to solving the relaxed problem is by using a direct method such as collocation or Bock's direct multiple shooting method. Optimal solutions for different control discretizations are also plotted in the leftmost figure.

The optimal objective value of this relaxed problem is <math>x_2(t_f) = 1.34408</math>. As follows from MIOC theory <bib id="Sager2011d" /> this is the best lower bound on the optimal value of the original problem with the integer restriction on the control function. In other words, this objective value can be approximated arbitrarily close, if the control only switches often enough between 0 and 1. As no optimal solution exists, two suboptimal ones are shown, one with only two switches and an objective function value of <math>x_2(t_f) = 1.38276</math>, and one with 56 switches and <math>x_2(t_f) = 1.34416</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:lotkaRelaxedControls.png| Optimal relaxed control determined by an indirect approach and by a direct approach on different control discretization grids.
Image:lotkaindirektStates.png| Differential states and corresponding adjoint variables in the indirect approach.
Image:lotka2Switches.png| Control and differential states with only two switches.
Image:lotka56Switches.png| Control and differential states with 56 switches.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:ACADO | ACADO code]] at [[Lotka Volterra fishing problem (ACADO)]]
* [[:Category:AMPL | AMPL code]] at [[Lotka Volterra fishing problem (AMPL)]]
* [[:Category:APMonitor | APMonitor code]] at [[Lotka Volterra fishing problem (APMonitor)]]
* [[:Category:Bocop | Bocop code]] at [[Lotka Volterra fishing problem (Bocop)]]
* [[:Category:Casadi | Casadi code]] at [[Lotka Volterra fishing problem (Casadi)]]
* [[:Category:Gekko | GEKKO Python code]] at [[Lotka Volterra fishing problem (GEKKO)]]
* [[:Category:JModelica | JModelica code]] at [[Lotka Volterra fishing problem (JModelica)]]
* [[:Category:Julia/JuMP | JuMP code]] at [[Lotka Volterra fishing problem (JuMP)]]
* [[:Category:Muscod | Muscod code]] at [[Lotka Volterra fishing problem (Muscod)]]
* [[:Category:switch | switch code]] at [[Lotka Volterra fishing problem (switch)]]
* [[:Category:TomDyn/PROPT | PROPT code]] at [[Lotka Volterra fishing problem (TomDyn/PROPT)]]

== Variants ==

There are several alternative formulations and variants of the above problem, in particular

* a prescribed time grid for the control function <bib id="Sager2006" />, see also [[Lotka Volterra fishing problem (AMPL)]],
* a time-optimal formulation to get into a steady-state <bib id="Sager2005" />,
* the usage of a different target steady-state, as the one corresponding to <math> w(t) = 1</math> which is <math>(1 + c_1, 1 - c_0)</math>, see [[Lotka Volterra multi-arcs problem]]
* different fishing control functions for the two species, see [[Lotka Volterra Multimode fishing problem]]
* different fishing control functions that fish an absolute value from the two species, see [[Lotka Volterra absolute fishing problem]]
* a terminal constrained formulation, where a violation is penalized via slack variables, see [[Lotka Volterra (terminal constraint violation)]]
* different parameters and start values.

== Miscellaneous and Further Reading ==
The Lotka Volterra fishing problem was introduced by Sebastian Sager in a proceedings paper <bib id="Sager2006" /> and revisited in his PhD thesis <bib id="Sager2005" />. These are also the references to look for more details.

== References ==
<biblist />


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

Quadrotor helicopter control problem

2019-10-14T13:11:06Z

ClemensZeile: Created page with "{{Dimensions |nd = 1 |nx = 6 |nu = 1 |nw = 3 |nre = 6 |nrc = 3 }}<!-- Do not insert line break here or Dimensions Box moves up in the l..."

{{Dimensions
|nd = 1
|nx = 6
|nu = 1
|nw = 3
|nre = 6
|nrc = 3
}}The mixed-integer optimal control problem of a quadrotor helicopter in two dimensions is taken from [7] and [8]. The evolution
of the quadrotor can be defined with respect to a fixed two dimensional reference
frame using six dimensions, where the first three dimensions represent the position
along a horizontal axis, the position along the vertical axis, and the roll angle of the
helicopter, respectively, and the last three dimensions represent the time derivative of
the first three dimensions.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w, s} & x_2(t_f) + 10s \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; c_0 x_0 \; w, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; c_1 x_1 \; w, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x_0 & \geq & 1.1 - s, \\
& x(0) &=& (0.5, 0.7, 0)^T, \\
& w(t) &\in& \{0, 1\}, \\
& s & \geq & 0.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. This problem variant penalizes a biomass x(0) that is below 1.1 at the end of the time horizon.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_{0}, c_{1}) &=& (0.4, 0.2),
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> is in the continuous interval <math>[0, 1]</math> rather than being binary, the optimal solution can be determined by means of direct optimal control.

The optimal objective value of the relaxed problem with <math> n_t=12000, \, n_u=200 </math> is <math>x_2(t_f) =1.36548113</math>. The objective value of the solution with binary controls obtained by Combinatorial Integral Approximation (CIA) is <math>x_2(t_f) =1.38756111</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:Lotka_term_ineq_Relaxed_12000_200.pdf| Optimal relaxed controls and states determined by an direct approach with ampl_mintoc (Radau collocation) and <math>n_t=12000, \, n_u=200</math>.
Image:Lotka_term_ineq_CIA_12000_200.pdf| Differential states and binary cotnrol determined by an direct approach (Radau collocation) with ampl_mintoc and <math>n_t=12000, \, n_u=200</math>. The relaxed controls were approximated by Combinatorial Integral Approximation.
</gallery>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

Lotka Volterra fishing problem

2019-10-14T13:07:26Z

ClemensZeile: /* Variants */

{{Dimensions
|nd = 1
|nx = 3
|nw = 1
|nre = 3
}}The '''Lotka Volterra fishing problem''' looks for an optimal fishing strategy to be performed on a fixed time horizon to bring the biomasses of both predator as prey fish to a prescribed steady state. The problem was set up as a small-scale benchmark problem.
The well known [http://en.wikipedia.org/wiki/Lotka_volterra Lotka Volterra equations] for a predator-prey system have been augmented by an additional linear term, relating to fishing by man. The control can be regarded both in a relaxed, as in a discrete manner, corresponding to a part of the fleet, or the full fishing fleet.

The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. The interior point equality conditions fix the initial values of the differential states.

The optimal integer control functions shows [[:Category:Chattering|chattering]] behavior, making the Lotka Volterra fishing problem an ideal candidate for benchmarking of algorithms.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; c_0 x_0 \; w, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; c_1 x_1 \; w, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x(0) &=& (0.5, 0.7, 0)^T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. The decision, whether the fishing fleet is actually fishing at time <math>t</math> is denoted by <math>w(t)</math>.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_0, c_1) &=& (0.4, 0.2).
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of [http://en.wikipedia.org/wiki/Pontryagin%27s_minimum_principle Pontryagins maximum principle]. The optimal solution contains a singular arc, as can be seen in the plot of the optimal control. The two differential states and corresponding adjoint variables in the indirect approach are also displayed. A different approach to solving the relaxed problem is by using a direct method such as collocation or Bock's direct multiple shooting method. Optimal solutions for different control discretizations are also plotted in the leftmost figure.

The optimal objective value of this relaxed problem is <math>x_2(t_f) = 1.34408</math>. As follows from MIOC theory <bib id="Sager2011d" /> this is the best lower bound on the optimal value of the original problem with the integer restriction on the control function. In other words, this objective value can be approximated arbitrarily close, if the control only switches often enough between 0 and 1. As no optimal solution exists, two suboptimal ones are shown, one with only two switches and an objective function value of <math>x_2(t_f) = 1.38276</math>, and one with 56 switches and <math>x_2(t_f) = 1.34416</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:lotkaRelaxedControls.png| Optimal relaxed control determined by an indirect approach and by a direct approach on different control discretization grids.
Image:lotkaindirektStates.png| Differential states and corresponding adjoint variables in the indirect approach.
Image:lotka2Switches.png| Control and differential states with only two switches.
Image:lotka56Switches.png| Control and differential states with 56 switches.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:ACADO | ACADO code]] at [[Lotka Volterra fishing problem (ACADO)]]
* [[:Category:AMPL | AMPL code]] at [[Lotka Volterra fishing problem (AMPL)]]
* [[:Category:APMonitor | APMonitor code]] at [[Lotka Volterra fishing problem (APMonitor)]]
* [[:Category:Bocop | Bocop code]] at [[Lotka Volterra fishing problem (Bocop)]]
* [[:Category:Casadi | Casadi code]] at [[Lotka Volterra fishing problem (Casadi)]]
* [[:Category:Gekko | GEKKO Python code]] at [[Lotka Volterra fishing problem (GEKKO)]]
* [[:Category:JModelica | JModelica code]] at [[Lotka Volterra fishing problem (JModelica)]]
* [[:Category:Julia/JuMP | JuMP code]] at [[Lotka Volterra fishing problem (JuMP)]]
* [[:Category:Muscod | Muscod code]] at [[Lotka Volterra fishing problem (Muscod)]]
* [[:Category:switch | switch code]] at [[Lotka Volterra fishing problem (switch)]]
* [[:Category:TomDyn/PROPT | PROPT code]] at [[Lotka Volterra fishing problem (TomDyn/PROPT)]]

== Variants ==

There are several alternative formulations and variants of the above problem, in particular

* a prescribed time grid for the control function <bib id="Sager2006" />, see also [[Lotka Volterra fishing problem (AMPL)]],
* a time-optimal formulation to get into a steady-state <bib id="Sager2005" />,
* the usage of a different target steady-state, as the one corresponding to <math> w(t) = 1</math> which is <math>(1 + c_1, 1 - c_0)</math>, see [[Lotka Volterra multi-arcs problem]]
* different fishing control functions for the two species, see [[Lotka Volterra Multimode fishing problem]]
* different fishing control functions that fish an absolute value from the two species, see [[Lotka Volterra absolute fishing problem]]
* a terminal constrained formulation, where a violation is penalized via slack variables, see [[Lotka Volterra (terminal constraint violation)]]
* different parameters and start values.
* [[Quadrotor helicopter control problem]]

== Miscellaneous and Further Reading ==
The Lotka Volterra fishing problem was introduced by Sebastian Sager in a proceedings paper <bib id="Sager2006" /> and revisited in his PhD thesis <bib id="Sager2005" />. These are also the references to look for more details.

== References ==
<biblist />


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

Lotka Volterra (terminal constraint violation)

2019-10-14T11:48:38Z

ClemensZeile: Created page with "{{Dimensions |nd = 1 |nx = 2 |nw = 1 |nre = 2 |nri = 1 }}Th..."

{{Dimensions
|nd = 1
|nx = 2
|nw = 1
|nre = 2
|nri = 1
}}This site describes a Lotka Volterra variant where a terminal inequality constraint on the differential states is added. A violation of this constraint is penalized as part of the Mayer objective.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w, s} & x_2(t_f) + 10s \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; c_0 x_0 \; w, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; c_1 x_1 \; w, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x_0 & \geq & 1.1 - s, \\
& x(0) &=& (0.5, 0.7, 0)^T, \\
& w(t) &\in& \{0, 1\}, \\
& s & \geq & 0.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. This problem variant penalizes a biomass x(0) that is below 1.1 at the end of the time horizon.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_{0}, c_{1}) &=& (0.4, 0.2),
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> is in the continuous interval <math>[0, 1]</math> rather than being binary, the optimal solution can be determined by means of direct optimal control.

The optimal objective value of the relaxed problem with <math> n_t=12000, \, n_u=200 </math> is <math>x_2(t_f) =1.36548113</math>. The objective value of the solution with binary controls obtained by Combinatorial Integral Approximation (CIA) is <math>x_2(t_f) =1.38756111</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:Lotka_term_ineq_Relaxed_12000_200.pdf| Optimal relaxed controls and states determined by an direct approach with ampl_mintoc (Radau collocation) and <math>n_t=12000, \, n_u=200</math>.
Image:Lotka_term_ineq_CIA_12000_200.pdf| Differential states and binary cotnrol determined by an direct approach (Radau collocation) with ampl_mintoc and <math>n_t=12000, \, n_u=200</math>. The relaxed controls were approximated by Combinatorial Integral Approximation.
</gallery>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

File:Lotka term ineq CIA 12000 200.pdf

2019-10-14T11:47:03Z

ClemensZeile: Numerical results for Lotka Volterra Problem with terminal inequality constraint violation. Results are based on binary control obtained via CIA and nt=12000, nu=200.

Numerical results for Lotka Volterra Problem with terminal inequality constraint violation. Results are based on binary control obtained via CIA and nt=12000, nu=200.

File:Lotka term ineq Relaxed 12000 200.pdf

2019-10-14T11:46:33Z

ClemensZeile: Numerical results for Lotka Volterra Problem with terminal inequality constraint violation. Results are based on relaxed control and nt=12000, nu=200.

Numerical results for Lotka Volterra Problem with terminal inequality constraint violation. Results are based on relaxed control and nt=12000, nu=200.

Lotka Volterra absolute fishing problem

2019-10-14T11:36:33Z

ClemensZeile: /* Mathematical formulation */

{{Dimensions
|nd = 1
|nx = 3
|nw = 5
|nre = 3
}}This site describes a Lotka Volterra variant with five binary controls that all represent fishing of an absolute biomass.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; \sum\limits_{i=1}^{5} c_{0,i}\; w_i, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; \sum\limits_{i=1}^{5} c_{1,i}\; w_i, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x(0) &=& (0.5, 0.7, 0)^T, \\
& \sum\limits_{i=1}^{5}w_i(t) &=& 1, \\
& w_i(t) &\in& \{0, 1\}, \quad i=1\ldots 5.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. This problem variant allows to choose between five different fishing options.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_{0,1}, c_{1,1}) &=& (0.2, 0.1),\\
(c_{0,2}, c_{1,2}) &=& (0.4, 0.2),\\
(c_{0,3}, c_{1,3}) &=& (0.01, 0.1),\\
(c_{0,4}, c_{1,4}) &=& (0, 0),\\
(c_{0,5}, c_{1,5}) &=& (-0.1, -0.2).
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> is in the continuous interval <math>[0, 1]</math> rather than being binary, the optimal solution can be determined by means of direct optimal control.

The optimal objective value of the relaxed problem with <math> n_t=12000, \, n_u=150 </math> is <math>x_2(t_f) =0.345768563</math>. The objective value of the solution with binary controls obtained by Combinatorial Integral Approximation (CIA) is <math>x_2(t_f) =0.348617982</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:Lotka_abs_fish_relaxed_12000_80.pdf| Optimal relaxed controls and states determined by an direct approach with ampl_mintoc (Radau collocation) and <math>n_t=12000, \, n_u=150</math>.
Image:Lotka_abs_fish_CIA_states_12000_80.pdf| Differential states determined by an direct approach (Radau collocation) with ampl_mintoc and <math>n_t=12000, \, n_u=150</math>. The relaxed controls were approximated by Combinatorial Integral Approximation.
Image:Lotka_abs_fish_CIA_controls_12000_80.pdf| Binary control determined by an direct approach (Radau collocation) with ampl_mintoc and <math>n_t=12000, \, n_u=150</math>. The relaxed controls were approximated by Combinatorial Integral Approximation. The fishing control shows a lot of chattering.
</gallery>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

Lotka Volterra fishing problem

2019-10-14T11:27:39Z

ClemensZeile: /* Variants */

{{Dimensions
|nd = 1
|nx = 3
|nw = 1
|nre = 3
}}The '''Lotka Volterra fishing problem''' looks for an optimal fishing strategy to be performed on a fixed time horizon to bring the biomasses of both predator as prey fish to a prescribed steady state. The problem was set up as a small-scale benchmark problem.
The well known [http://en.wikipedia.org/wiki/Lotka_volterra Lotka Volterra equations] for a predator-prey system have been augmented by an additional linear term, relating to fishing by man. The control can be regarded both in a relaxed, as in a discrete manner, corresponding to a part of the fleet, or the full fishing fleet.

The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. The interior point equality conditions fix the initial values of the differential states.

The optimal integer control functions shows [[:Category:Chattering|chattering]] behavior, making the Lotka Volterra fishing problem an ideal candidate for benchmarking of algorithms.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; c_0 x_0 \; w, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; c_1 x_1 \; w, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x(0) &=& (0.5, 0.7, 0)^T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. The decision, whether the fishing fleet is actually fishing at time <math>t</math> is denoted by <math>w(t)</math>.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_0, c_1) &=& (0.4, 0.2).
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of [http://en.wikipedia.org/wiki/Pontryagin%27s_minimum_principle Pontryagins maximum principle]. The optimal solution contains a singular arc, as can be seen in the plot of the optimal control. The two differential states and corresponding adjoint variables in the indirect approach are also displayed. A different approach to solving the relaxed problem is by using a direct method such as collocation or Bock's direct multiple shooting method. Optimal solutions for different control discretizations are also plotted in the leftmost figure.

The optimal objective value of this relaxed problem is <math>x_2(t_f) = 1.34408</math>. As follows from MIOC theory <bib id="Sager2011d" /> this is the best lower bound on the optimal value of the original problem with the integer restriction on the control function. In other words, this objective value can be approximated arbitrarily close, if the control only switches often enough between 0 and 1. As no optimal solution exists, two suboptimal ones are shown, one with only two switches and an objective function value of <math>x_2(t_f) = 1.38276</math>, and one with 56 switches and <math>x_2(t_f) = 1.34416</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:lotkaRelaxedControls.png| Optimal relaxed control determined by an indirect approach and by a direct approach on different control discretization grids.
Image:lotkaindirektStates.png| Differential states and corresponding adjoint variables in the indirect approach.
Image:lotka2Switches.png| Control and differential states with only two switches.
Image:lotka56Switches.png| Control and differential states with 56 switches.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:ACADO | ACADO code]] at [[Lotka Volterra fishing problem (ACADO)]]
* [[:Category:AMPL | AMPL code]] at [[Lotka Volterra fishing problem (AMPL)]]
* [[:Category:APMonitor | APMonitor code]] at [[Lotka Volterra fishing problem (APMonitor)]]
* [[:Category:Bocop | Bocop code]] at [[Lotka Volterra fishing problem (Bocop)]]
* [[:Category:Casadi | Casadi code]] at [[Lotka Volterra fishing problem (Casadi)]]
* [[:Category:Gekko | GEKKO Python code]] at [[Lotka Volterra fishing problem (GEKKO)]]
* [[:Category:JModelica | JModelica code]] at [[Lotka Volterra fishing problem (JModelica)]]
* [[:Category:Julia/JuMP | JuMP code]] at [[Lotka Volterra fishing problem (JuMP)]]
* [[:Category:Muscod | Muscod code]] at [[Lotka Volterra fishing problem (Muscod)]]
* [[:Category:switch | switch code]] at [[Lotka Volterra fishing problem (switch)]]
* [[:Category:TomDyn/PROPT | PROPT code]] at [[Lotka Volterra fishing problem (TomDyn/PROPT)]]

== Variants ==

There are several alternative formulations and variants of the above problem, in particular

* a prescribed time grid for the control function <bib id="Sager2006" />, see also [[Lotka Volterra fishing problem (AMPL)]],
* a time-optimal formulation to get into a steady-state <bib id="Sager2005" />,
* the usage of a different target steady-state, as the one corresponding to <math> w(t) = 1</math> which is <math>(1 + c_1, 1 - c_0)</math>, see [[Lotka Volterra multi-arcs problem]]
* different fishing control functions for the two species, see [[Lotka Volterra Multimode fishing problem]]
* different fishing control functions that fish an absolute value from the two species, see [[Lotka Volterra absolute fishing problem]]
* a terminal constrained formulation, where a violation is penalized via slack variables, see [[Lotka Volterra (terminal constraint violation)]]
* different parameters and start values.

== Miscellaneous and Further Reading ==
The Lotka Volterra fishing problem was introduced by Sebastian Sager in a proceedings paper <bib id="Sager2006" /> and revisited in his PhD thesis <bib id="Sager2005" />. These are also the references to look for more details.

== References ==
<biblist />


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

Lotka Volterra absolute fishing problem

2019-10-14T11:10:19Z

ClemensZeile: /* Reference Solutions */

{{Dimensions
|nd = 1
|nx = 3
|nw = 5
|nre = 3
}}This site describes a Lotka Volterra variant with five binary controls that all represent fishing of an absolute biomass.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; \sum\limits_{i=1}^{5} c_{0,i}\; w_i, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; \sum\limits_{i=1}^{5} c_{1,i}\; w_i, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x(0) &=& (0.5, 0.7, 0)^T, \\
& \sum\limits_{i=1}^{5}w_i(t) &=& 1, \\
& w_i(t) &\in& \{0, 1\}, \quad i=1\ldots 5.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. This problem variant allows to choose between three different fishing options.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_{0,1}, c_{1,1}) &=& (0.2, 0.1),\\
(c_{0,2}, c_{1,2}) &=& (0.4, 0.2),\\
(c_{0,3}, c_{1,3}) &=& (0.01, 0.1),\\
(c_{0,4}, c_{1,4}) &=& (0, 0),\\
(c_{0,5}, c_{1,5}) &=& (-0.1, -0.2).
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> is in the continuous interval <math>[0, 1]</math> rather than being binary, the optimal solution can be determined by means of direct optimal control.

The optimal objective value of the relaxed problem with <math> n_t=12000, \, n_u=150 </math> is <math>x_2(t_f) =0.345768563</math>. The objective value of the solution with binary controls obtained by Combinatorial Integral Approximation (CIA) is <math>x_2(t_f) =0.348617982</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:Lotka_abs_fish_relaxed_12000_80.pdf| Optimal relaxed controls and states determined by an direct approach with ampl_mintoc (Radau collocation) and <math>n_t=12000, \, n_u=150</math>.
Image:Lotka_abs_fish_CIA_states_12000_80.pdf| Differential states determined by an direct approach (Radau collocation) with ampl_mintoc and <math>n_t=12000, \, n_u=150</math>. The relaxed controls were approximated by Combinatorial Integral Approximation.
Image:Lotka_abs_fish_CIA_controls_12000_80.pdf| Binary control determined by an direct approach (Radau collocation) with ampl_mintoc and <math>n_t=12000, \, n_u=150</math>. The relaxed controls were approximated by Combinatorial Integral Approximation. The fishing control shows a lot of chattering.
</gallery>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

File:Lotka abs fish CIA controls 12000 80.pdf

2019-10-14T11:08:34Z

ClemensZeile: Results of Lotka Volterra Problem with absolute fishing controls obtained via CIA, control values

Results of Lotka Volterra Problem with absolute fishing controls obtained via CIA, control values

File:Lotka abs fish CIA states 12000 80.pdf

2019-10-14T11:08:14Z

ClemensZeile: Results of Lotka Volterra Problem with absolute fishing controls obtained via CIA, differential states

Results of Lotka Volterra Problem with absolute fishing controls obtained via CIA, differential states

File:Lotka abs fish relaxed 12000 80.pdf

2019-10-14T11:03:51Z

ClemensZeile: Results of Lotka Volterra absolute fishing problem with relaxed controls.

Results of Lotka Volterra absolute fishing problem with relaxed controls.

Lotka Volterra absolute fishing problem

2019-10-14T10:56:25Z

ClemensZeile:

{{Dimensions
|nd = 1
|nx = 3
|nw = 5
|nre = 3
}}This site describes a Lotka Volterra variant with five binary controls that all represent fishing of an absolute biomass.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; \sum\limits_{i=1}^{5} c_{0,i}\; w_i, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; \sum\limits_{i=1}^{5} c_{1,i}\; w_i, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x(0) &=& (0.5, 0.7, 0)^T, \\
& \sum\limits_{i=1}^{5}w_i(t) &=& 1, \\
& w_i(t) &\in& \{0, 1\}, \quad i=1\ldots 5.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. This problem variant allows to choose between three different fishing options.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_{0,1}, c_{1,1}) &=& (0.2, 0.1),\\
(c_{0,2}, c_{1,2}) &=& (0.4, 0.2),\\
(c_{0,3}, c_{1,3}) &=& (0.01, 0.1),\\
(c_{0,4}, c_{1,4}) &=& (0, 0),\\
(c_{0,5}, c_{1,5}) &=& (-0.1, -0.2).
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> is in the continuous interval <math>[0, 1]</math> rather than being binary, the optimal solution can be determined by means of direct optimal control.

The optimal objective value of the relaxed problem with <math> n_t=12000, \, n_u=150 </math> is <math>x_2(t_f) =0.345768563</math>. The objective value of the solution with binary controls obtained by Combinatorial Integral Approximation (CIA) is <math>x_2(t_f) =0.348617982</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:MmlotkaRelaxed_12000_30_1.png| Optimal relaxed controls and states determined by an direct approach with ampl_mintoc (Radau collocation) and <math>n_t=12000, \, n_u=400</math>.
Image:MmlotkaCIA 12000 30 1.png| Optimal binary controls and states determined by an direct approach (Radau collocation) with ampl_mintoc and <math>n_t=12000, \, n_u=400</math>. The relaxed controls were approximated by Combinatorial Integral Approximation.
</gallery>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

Lotka Volterra absolute fishing problem

2019-10-14T07:24:38Z

ClemensZeile: Created page with "{{Dimensions |nd = 1 |nx = 3 |nw = 5 |nre = 3 }}This site describ..."

{{Dimensions
|nd = 1
|nx = 3
|nw = 5
|nre = 3
}}This site describes a Lotka Volterra variant with five binary controls that all represent fishing of an absolute biomass.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; \sum\limits_{i=1}^{5} c_{0,i}\; w_i, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; \sum\limits_{i=1}^{5} c_{1,i}\; w_i, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x(0) &=& (0.5, 0.7, 0)^T, \\
& \sum\limits_{i=1}^{5}w_i(t) &=& 1, \\
& w_i(t) &\in& \{0, 1\}, \quad i=1\ldots 5.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. This problem variant allows to choose between three different fishing options.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_{0,1}, c_{1,1}) &=& (0.2, 0.1),\\
(c_{0,2}, c_{1,2}) &=& (0.4, 0.2),\\
(c_{0,3}, c_{1,3}) &=& (0.01, 0.1),\\
(c_{0,4}, c_{1,4}) &=& (0, 0),\\
(c_{0,5}, c_{1,5}) &=& (-0.1, -0.2).
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of direct optimal control.

The optimal objective value of the relaxed problem with <math> n_t=12000, \, n_u=400 </math> is <math>x_2(t_f) =1.82875272</math>. The objective value of the binary controls obtained by Combinatorial Integral Approimation (CIA) is <math>x_2(t_f) =1.82878681</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:MmlotkaRelaxed_12000_30_1.png| Optimal relaxed controls and states determined by an direct approach with ampl_mintoc (Radau collocation) and <math>n_t=12000, \, n_u=400</math>.
Image:MmlotkaCIA 12000 30 1.png| Optimal binary controls and states determined by an direct approach (Radau collocation) with ampl_mintoc and <math>n_t=12000, \, n_u=400</math>. The relaxed controls were approximated by Combinatorial Integral Approximation.
</gallery>


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

Lotka Volterra fishing problem

2019-10-14T06:34:26Z

ClemensZeile: /* Variants */

{{Dimensions
|nd = 1
|nx = 3
|nw = 1
|nre = 3
}}The '''Lotka Volterra fishing problem''' looks for an optimal fishing strategy to be performed on a fixed time horizon to bring the biomasses of both predator as prey fish to a prescribed steady state. The problem was set up as a small-scale benchmark problem.
The well known [http://en.wikipedia.org/wiki/Lotka_volterra Lotka Volterra equations] for a predator-prey system have been augmented by an additional linear term, relating to fishing by man. The control can be regarded both in a relaxed, as in a discrete manner, corresponding to a part of the fleet, or the full fishing fleet.

The mathematical equations form a small-scale [[:Category:ODE model|ODE model]]. The interior point equality conditions fix the initial values of the differential states.

The optimal integer control functions shows [[:Category:Chattering|chattering]] behavior, making the Lotka Volterra fishing problem an ideal candidate for benchmarking of algorithms.

== Mathematical formulation ==

The mixed-integer optimal control problem is given by


<math>
\begin{array}{llclr}
\displaystyle \min_{x, w} & x_2(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_0 & = & x_0 - x_0 x_1 - \; c_0 x_0 \; w, \\
& \dot{x}_1 & = & - x_1 + x_0 x_1 - \; c_1 x_1 \; w, \\
& \dot{x}_2 & = & (x_0 - 1)^2 + (x_1 - 1)^2, \\[1.5ex]
& x(0) &=& (0.5, 0.7, 0)^T, \\
& w(t) &\in& \{0, 1\}.
\end{array}
</math>


Here the differential states <math>(x_0, x_1)</math> describe the biomasses of prey and predator, respectively. The third differential state is used here to transform the objective, an integrated deviation, into the Mayer formulation <math>\min \; x_2(t_f)</math>. The decision, whether the fishing fleet is actually fishing at time <math>t</math> is denoted by <math>w(t)</math>.

== Parameters ==

These fixed values are used within the model.

<math>
\begin{array}{rcl}
[t_0, t_f] &=& [0, 12],\\
(c_0, c_1) &=& (0.4, 0.2).
\end{array}
</math>

== Reference Solutions ==

If the problem is relaxed, i.e., we demand that <math>w(t)</math> be in the continuous interval <math>[0, 1]</math> instead of the binary choice <math>\{0,1\}</math>, the optimal solution can be determined by means of [http://en.wikipedia.org/wiki/Pontryagin%27s_minimum_principle Pontryagins maximum principle]. The optimal solution contains a singular arc, as can be seen in the plot of the optimal control. The two differential states and corresponding adjoint variables in the indirect approach are also displayed. A different approach to solving the relaxed problem is by using a direct method such as collocation or Bock's direct multiple shooting method. Optimal solutions for different control discretizations are also plotted in the leftmost figure.

The optimal objective value of this relaxed problem is <math>x_2(t_f) = 1.34408</math>. As follows from MIOC theory <bib id="Sager2011d" /> this is the best lower bound on the optimal value of the original problem with the integer restriction on the control function. In other words, this objective value can be approximated arbitrarily close, if the control only switches often enough between 0 and 1. As no optimal solution exists, two suboptimal ones are shown, one with only two switches and an objective function value of <math>x_2(t_f) = 1.38276</math>, and one with 56 switches and <math>x_2(t_f) = 1.34416</math>.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:lotkaRelaxedControls.png| Optimal relaxed control determined by an indirect approach and by a direct approach on different control discretization grids.
Image:lotkaindirektStates.png| Differential states and corresponding adjoint variables in the indirect approach.
Image:lotka2Switches.png| Control and differential states with only two switches.
Image:lotka56Switches.png| Control and differential states with 56 switches.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:ACADO | ACADO code]] at [[Lotka Volterra fishing problem (ACADO)]]
* [[:Category:AMPL | AMPL code]] at [[Lotka Volterra fishing problem (AMPL)]]
* [[:Category:APMonitor | APMonitor code]] at [[Lotka Volterra fishing problem (APMonitor)]]
* [[:Category:Bocop | Bocop code]] at [[Lotka Volterra fishing problem (Bocop)]]
* [[:Category:Casadi | Casadi code]] at [[Lotka Volterra fishing problem (Casadi)]]
* [[:Category:Gekko | GEKKO Python code]] at [[Lotka Volterra fishing problem (GEKKO)]]
* [[:Category:JModelica | JModelica code]] at [[Lotka Volterra fishing problem (JModelica)]]
* [[:Category:Julia/JuMP | JuMP code]] at [[Lotka Volterra fishing problem (JuMP)]]
* [[:Category:Muscod | Muscod code]] at [[Lotka Volterra fishing problem (Muscod)]]
* [[:Category:switch | switch code]] at [[Lotka Volterra fishing problem (switch)]]
* [[:Category:TomDyn/PROPT | PROPT code]] at [[Lotka Volterra fishing problem (TomDyn/PROPT)]]

== Variants ==

There are several alternative formulations and variants of the above problem, in particular

* a prescribed time grid for the control function <bib id="Sager2006" />, see also [[Lotka Volterra fishing problem (AMPL)]],
* a time-optimal formulation to get into a steady-state <bib id="Sager2005" />,
* the usage of a different target steady-state, as the one corresponding to <math> w(t) = 1</math> which is <math>(1 + c_1, 1 - c_0)</math>, see [[Lotka Volterra multi-arcs problem]]
* different fishing control functions for the two species, see [[Lotka Volterra Multimode fishing problem]]
* different fishing control functions that fish an absolute value from the two species, see [[Lotka Volterra absolute fishing problem]]
* different parameters and start values.

== Miscellaneous and Further Reading ==
The Lotka Volterra fishing problem was introduced by Sebastian Sager in a proceedings paper <bib id="Sager2006" /> and revisited in his PhD thesis <bib id="Sager2005" />. These are also the references to look for more details.

== References ==
<biblist />


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Tracking objective]]
[[Category:Chattering]]
[[Category:Sensitivity-seeking arcs]]
[[Category:Population dynamics]]

Continuously Stirred Tank Reactor problem

2019-09-02T10:11:39Z

ClemensZeile: /* Parameters */

{{Dimensions
|nd = 1
|nx = 4
|nu = 2
|nre = 2
}}The Continuously Stirred Tank Reactor problem considers a chemical reaction that produces cyclopenthenol while using up cyclepentadiene "by an acid-catalyzed electrophilic hydration in aqueous solution", an exothermal reaction which needs to be cooled. This problem can e.g. be found in <bib id="Diehl2001" />.

The inflow into the tank contains only cyclopentadiene (substance <math> A </math>) with temperature <math> \theta_0 </math> and the flow rate <math> \dot{V} </math> can be controlled. The outflow rate is the same as the inflow rate to keep the liquid level in the tank constant.
"The outflow contains a remainder of cyclopentadiene, the wanted product cyclepentenol (substance <math> B </math>) and two unwated by-products, cyclopentanediol (substance <math> C </math>) and dicyclopentadiene (substance <math> D </math>) with concentrations <math> c_A, c_B, c_C, c_D </math>." The latter two are not tracked in the problem as the substances are not of use.
The reaction scheme is given as:

<math>
\begin{array}{ccccc}
A & \overset{k_1}\rightarrow & B \overset{k_2}\rightarrow & C\\
2A & \overset{k_3}\rightarrow & D
\end{array}
</math>

where the reaction rates <math> k_i </math> are a function of the reactor temperature <math> \theta </math> via an Arrhenius law
<math> k_i(\theta) = k_{i0} \cdot \exp ( \frac{E_i}{\theta / ^\circ C + 273.15} ), \quad i=1,2,3. </math>

"The temperature <math> \theta_K </math> in the cooling jacket is held down by an external heat exchanger whose heat removal rate <math> \dot{Q}_K </math> can be controlled."

== Mathematical formulation ==
The problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{\dot{V}, \dot{Q}_K} & c_B & & \text{ at the end of reaction} \\[1.5ex]
\mbox{s.t.}
& \dot{c_A} & = & \frac{\dot{V}}{V_R} (c_{A0} - c_A) - k_1 c_A - k_3 c_A^2, \\[0.5cm]
& \dot{c_B} & = & -\frac{\dot{V}}{V_R} c_B + k_1 c_A - k_2 c_B, \\[0.5cm]
& \dot{\theta} & = & \frac{\dot{V}}{V_R} ( \theta_0 - \theta) + \frac{k_w A_R}{\rho C_p V_R} (\theta_K - \theta) - \frac{1}{\rho C_p} (k_1 c_A H_1 + k_2 c_B H_2 + k_3 c_A^2 H_3), \\[0.5cm]
& \dot{\theta_K} & = & \frac{1}{m_K C_{PK}} ( \dot{Q}_K + k_w A_R (\theta - \theta_K)),\\[0.7cm]
& c_A(0) & = & c_{A0},\\
& c_B(0) & = & 0.
\end{array}
</math>

where the various values are given in the Parameters section.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Arrhenius coefficient
|<math>k_{10}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{20}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{30}</math>
|<math> 9.043 \cdot 10^{9} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>E_1</math>
|<math> -9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_2</math>
|<math>-9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_3</math>
|<math>-8560</math>
|[-]
|-
|Reaction enthalpy
|<math>H_1</math>
|<math>4.2</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_2</math>
|<math>-11.0</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_3</math>
|<math>-41.85</math>
|<math> \frac{kJ}{mol} </math>
|-
|Solution density
|<math>\rho</math>
|<math>0.9342</math>
|<math> \frac{kg}{l} </math>
|-
|Capacity of aqueous solution
|<math>C_p</math>
|<math>3.01</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Heat transfer coefficient for cooling jacket
|<math>k_w</math>
|<math>4032</math>
|<math> \frac{kJ}{h \cdot m^2 \cdot K} </math>
|-
|Reactor surface area
|<math>A_R</math>
|<math>0.215</math>
|<math> m^2 </math>
|-
|Reactor volume
|<math>V_R</math>
|<math>10</math>
|<math> l </math>
|-
|Coolant mass
|<math>m_K</math>
|<math>5</math>
|<math> kg </math>
|-
|Capacity of coolant solution
|<math>C_{PK}</math>
|<math>2.0</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Starting concentration of subs. <math> A </math>
|<math>c_{A0}</math>
|<math>5.1</math>
|<math> \frac{mol}{l} </math>
|-
|Inflow temperature
|<math>\theta_0</math>
|<math>104.9</math>
|<math> ^\circ C </math>
|}

== Reference solution ==

"The result of a steady state optimization of the yield <math> = \frac{c_B |_S}{c_{A0}} </math> with respect to the design parameter <math> \theta_0 </math> (feed temperature) and the two controls yields the steady stae and controls"
<math> c_A = 2.1402 \frac{mol}{l}, c_B = 1.0903\frac{mol}{l}, \theta = 114.19^\circ C, \theta_K = 112.91^\circ C </math> and <math> \frac{\dot{V}}{V_R} = 14.19 h^{-1}, \dot{Q}_K = -1113.5 \frac{kJ}{h} </math>.

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Continuously Stirred Tank Reactor (TACO)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Continuously Stirred Tank Reactor problem

2019-09-02T10:11:16Z

ClemensZeile: /* Parameters */

{{Dimensions
|nd = 1
|nx = 4
|nu = 2
|nre = 2
}}The Continuously Stirred Tank Reactor problem considers a chemical reaction that produces cyclopenthenol while using up cyclepentadiene "by an acid-catalyzed electrophilic hydration in aqueous solution", an exothermal reaction which needs to be cooled. This problem can e.g. be found in <bib id="Diehl2001" />.

The inflow into the tank contains only cyclopentadiene (substance <math> A </math>) with temperature <math> \theta_0 </math> and the flow rate <math> \dot{V} </math> can be controlled. The outflow rate is the same as the inflow rate to keep the liquid level in the tank constant.
"The outflow contains a remainder of cyclopentadiene, the wanted product cyclepentenol (substance <math> B </math>) and two unwated by-products, cyclopentanediol (substance <math> C </math>) and dicyclopentadiene (substance <math> D </math>) with concentrations <math> c_A, c_B, c_C, c_D </math>." The latter two are not tracked in the problem as the substances are not of use.
The reaction scheme is given as:

<math>
\begin{array}{ccccc}
A & \overset{k_1}\rightarrow & B \overset{k_2}\rightarrow & C\\
2A & \overset{k_3}\rightarrow & D
\end{array}
</math>

where the reaction rates <math> k_i </math> are a function of the reactor temperature <math> \theta </math> via an Arrhenius law
<math> k_i(\theta) = k_{i0} \cdot \exp ( \frac{E_i}{\theta / ^\circ C + 273.15} ), \quad i=1,2,3. </math>

"The temperature <math> \theta_K </math> in the cooling jacket is held down by an external heat exchanger whose heat removal rate <math> \dot{Q}_K </math> can be controlled."

== Mathematical formulation ==
The problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{\dot{V}, \dot{Q}_K} & c_B & & \text{ at the end of reaction} \\[1.5ex]
\mbox{s.t.}
& \dot{c_A} & = & \frac{\dot{V}}{V_R} (c_{A0} - c_A) - k_1 c_A - k_3 c_A^2, \\[0.5cm]
& \dot{c_B} & = & -\frac{\dot{V}}{V_R} c_B + k_1 c_A - k_2 c_B, \\[0.5cm]
& \dot{\theta} & = & \frac{\dot{V}}{V_R} ( \theta_0 - \theta) + \frac{k_w A_R}{\rho C_p V_R} (\theta_K - \theta) - \frac{1}{\rho C_p} (k_1 c_A H_1 + k_2 c_B H_2 + k_3 c_A^2 H_3), \\[0.5cm]
& \dot{\theta_K} & = & \frac{1}{m_K C_{PK}} ( \dot{Q}_K + k_w A_R (\theta - \theta_K)),\\[0.7cm]
& c_A(0) & = & c_{A0},\\
& c_B(0) & = & 0.
\end{array}
</math>

where the various values are given in the Parameters section.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Arrhenius coefficient
|<math>k_{10}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{20}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{30}</math>
|<math> 9.043 \cdot 10 </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>E_1</math>
|<math> -9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_2</math>
|<math>-9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_3</math>
|<math>-8560</math>
|[-]
|-
|Reaction enthalpy
|<math>H_1</math>
|<math>4.2</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_2</math>
|<math>-11.0</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_3</math>
|<math>-41.85</math>
|<math> \frac{kJ}{mol} </math>
|-
|Solution density
|<math>\rho</math>
|<math>0.9342</math>
|<math> \frac{kg}{l} </math>
|-
|Capacity of aqueous solution
|<math>C_p</math>
|<math>3.01</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Heat transfer coefficient for cooling jacket
|<math>k_w</math>
|<math>4032</math>
|<math> \frac{kJ}{h \cdot m^2 \cdot K} </math>
|-
|Reactor surface area
|<math>A_R</math>
|<math>0.215</math>
|<math> m^2 </math>
|-
|Reactor volume
|<math>V_R</math>
|<math>10</math>
|<math> l </math>
|-
|Coolant mass
|<math>m_K</math>
|<math>5</math>
|<math> kg </math>
|-
|Capacity of coolant solution
|<math>C_{PK}</math>
|<math>2.0</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Starting concentration of subs. <math> A </math>
|<math>c_{A0}</math>
|<math>5.1</math>
|<math> \frac{mol}{l} </math>
|-
|Inflow temperature
|<math>\theta_0</math>
|<math>104.9</math>
|<math> ^\circ C </math>
|}

== Reference solution ==

"The result of a steady state optimization of the yield <math> = \frac{c_B |_S}{c_{A0}} </math> with respect to the design parameter <math> \theta_0 </math> (feed temperature) and the two controls yields the steady stae and controls"
<math> c_A = 2.1402 \frac{mol}{l}, c_B = 1.0903\frac{mol}{l}, \theta = 114.19^\circ C, \theta_K = 112.91^\circ C </math> and <math> \frac{\dot{V}}{V_R} = 14.19 h^{-1}, \dot{Q}_K = -1113.5 \frac{kJ}{h} </math>.

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Continuously Stirred Tank Reactor (TACO)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Continuously Stirred Tank Reactor problem

2019-09-02T10:02:31Z

ClemensZeile: /* Parameters */

Continuously Stirred Tank Reactor problem

2019-09-02T09:52:45Z

ClemensZeile: /* Parameters */

{{Dimensions
|nd = 1
|nx = 4
|nu = 2
|nre = 2
}}The Continuously Stirred Tank Reactor problem considers a chemical reaction that produces cyclopenthenol while using up cyclepentadiene "by an acid-catalyzed electrophilic hydration in aqueous solution", an exothermal reaction which needs to be cooled. This problem can e.g. be found in <bib id="Diehl2001" />.

The inflow into the tank contains only cyclopentadiene (substance <math> A </math>) with temperature <math> \theta_0 </math> and the flow rate <math> \dot{V} </math> can be controlled. The outflow rate is the same as the inflow rate to keep the liquid level in the tank constant.
"The outflow contains a remainder of cyclopentadiene, the wanted product cyclepentenol (substance <math> B </math>) and two unwated by-products, cyclopentanediol (substance <math> C </math>) and dicyclopentadiene (substance <math> D </math>) with concentrations <math> c_A, c_B, c_C, c_D </math>." The latter two are not tracked in the problem as the substances are not of use.
The reaction scheme is given as:

<math>
\begin{array}{ccccc}
A & \overset{k_1}\rightarrow & B \overset{k_2}\rightarrow & C\\
2A & \overset{k_3}\rightarrow & D
\end{array}
</math>

where the reaction rates <math> k_i </math> are a function of the reactor temperature <math> \theta </math> via an Arrhenius law
<math> k_i(\theta) = k_{i0} \cdot \exp ( \frac{E_i}{\theta / ^\circ C + 273.15} ), \quad i=1,2,3. </math>

"The temperature <math> \theta_K </math> in the cooling jacket is held down by an external heat exchanger whose heat removal rate <math> \dot{Q}_K </math> can be controlled."

== Mathematical formulation ==
The problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{\dot{V}, \dot{Q}_K} & c_B & & \text{ at the end of reaction} \\[1.5ex]
\mbox{s.t.}
& \dot{c_A} & = & \frac{\dot{V}}{V_R} (c_{A0} - c_A) - k_1 c_A - k_3 c_A^2, \\[0.5cm]
& \dot{c_B} & = & -\frac{\dot{V}}{V_R} c_B + k_1 c_A - k_2 c_B, \\[0.5cm]
& \dot{\theta} & = & \frac{\dot{V}}{V_R} ( \theta_0 - \theta) + \frac{k_w A_R}{\rho C_p V_R} (\theta_K - \theta) - \frac{1}{\rho C_p} (k_1 c_A H_1 + k_2 c_B H_2 + k_3 c_A^2 H_3), \\[0.5cm]
& \dot{\theta_K} & = & \frac{1}{m_K C_{PK}} ( \dot{Q}_K + k_w A_R (\theta - \theta_K)),\\[0.7cm]
& c_A(0) & = & c_{A0},\\
& c_B(0) & = & 0.
\end{array}
</math>

where the various values are given in the Parameters section.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Arrhenius coefficient
|<math>k_{10}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{20}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{30}</math>
|<math> 9.043 \cdot 10^{19} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>E_1</math>
|<math> -9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_2</math>
|<math>-9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_3</math>
|<math>-8560</math>
|[-]
|-
|Reaction enthalpy
|<math>H_1</math>
|<math>4.2</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_2</math>
|<math>-11.0</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_3</math>
|<math>-41.85</math>
|<math> \frac{kJ}{mol} </math>
|-
|Solution density
|<math>\rho</math>
|<math>0.9342</math>
|<math> \frac{kg}{l} </math>
|-
|Capacity of aqueous solution
|<math>C_p</math>
|<math>3.01</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Heat transfer coefficient for cooling jacket
|<math>k_w</math>
|<math>4032</math>
|<math> \frac{kJ}{h \cdot m^2 \cdot K} </math>
|-
|Reactor surface area
|<math>A_R</math>
|<math>0.215</math>
|<math> m^2 </math>
|-
|Reactor volume
|<math>V_R</math>
|<math>10</math>
|<math> l </math>
|-
|Coolant mass
|<math>m_K</math>
|<math>5</math>
|<math> kg </math>
|-
|Capacity of coolant solution
|<math>C_{PK}</math>
|<math>2.0</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Starting concentration of subs. <math> A </math>
|<math>c_{A0}</math>
|<math>5.1</math>
|<math> \frac{mol}{l} </math>
|-
|Inflow temperature
|<math>\theta_0</math>
|<math>104.9</math>
|<math> ^\circ C </math>
|}

== Reference solution ==

"The result of a steady state optimization of the yield <math> = \frac{c_B |_S}{c_{A0}} </math> with respect to the design parameter <math> \theta_0 </math> (feed temperature) and the two controls yields the steady stae and controls"
<math> c_A = 2.1402 \frac{mol}{l}, c_B = 1.0903\frac{mol}{l}, \theta = 114.19^\circ C, \theta_K = 112.91^\circ C </math> and <math> \frac{\dot{V}}{V_R} = 14.19 h^{-1}, \dot{Q}_K = -1113.5 \frac{kJ}{h} </math>.

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Continuously Stirred Tank Reactor (TACO)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Continuously Stirred Tank Reactor problem

2019-09-02T09:52:18Z

ClemensZeile: /* Parameters */

{{Dimensions
|nd = 1
|nx = 4
|nu = 2
|nre = 2
}}The Continuously Stirred Tank Reactor problem considers a chemical reaction that produces cyclopenthenol while using up cyclepentadiene "by an acid-catalyzed electrophilic hydration in aqueous solution", an exothermal reaction which needs to be cooled. This problem can e.g. be found in <bib id="Diehl2001" />.

The inflow into the tank contains only cyclopentadiene (substance <math> A </math>) with temperature <math> \theta_0 </math> and the flow rate <math> \dot{V} </math> can be controlled. The outflow rate is the same as the inflow rate to keep the liquid level in the tank constant.
"The outflow contains a remainder of cyclopentadiene, the wanted product cyclepentenol (substance <math> B </math>) and two unwated by-products, cyclopentanediol (substance <math> C </math>) and dicyclopentadiene (substance <math> D </math>) with concentrations <math> c_A, c_B, c_C, c_D </math>." The latter two are not tracked in the problem as the substances are not of use.
The reaction scheme is given as:

<math>
\begin{array}{ccccc}
A & \overset{k_1}\rightarrow & B \overset{k_2}\rightarrow & C\\
2A & \overset{k_3}\rightarrow & D
\end{array}
</math>

where the reaction rates <math> k_i </math> are a function of the reactor temperature <math> \theta </math> via an Arrhenius law
<math> k_i(\theta) = k_{i0} \cdot \exp ( \frac{E_i}{\theta / ^\circ C + 273.15} ), \quad i=1,2,3. </math>

"The temperature <math> \theta_K </math> in the cooling jacket is held down by an external heat exchanger whose heat removal rate <math> \dot{Q}_K </math> can be controlled."

== Mathematical formulation ==
The problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{\dot{V}, \dot{Q}_K} & c_B & & \text{ at the end of reaction} \\[1.5ex]
\mbox{s.t.}
& \dot{c_A} & = & \frac{\dot{V}}{V_R} (c_{A0} - c_A) - k_1 c_A - k_3 c_A^2, \\[0.5cm]
& \dot{c_B} & = & -\frac{\dot{V}}{V_R} c_B + k_1 c_A - k_2 c_B, \\[0.5cm]
& \dot{\theta} & = & \frac{\dot{V}}{V_R} ( \theta_0 - \theta) + \frac{k_w A_R}{\rho C_p V_R} (\theta_K - \theta) - \frac{1}{\rho C_p} (k_1 c_A H_1 + k_2 c_B H_2 + k_3 c_A^2 H_3), \\[0.5cm]
& \dot{\theta_K} & = & \frac{1}{m_K C_{PK}} ( \dot{Q}_K + k_w A_R (\theta - \theta_K)),\\[0.7cm]
& c_A(0) & = & c_{A0},\\
& c_B(0) & = & 0.
\end{array}
</math>

where the various values are given in the Parameters section.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Arrhenius coefficient
|<math>k_{10}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{20}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{30}</math>
|<math> 9.043 \cdot 10^{9\ } </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>E_1</math>
|<math> -9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_2</math>
|<math>-9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_3</math>
|<math>-8560</math>
|[-]
|-
|Reaction enthalpy
|<math>H_1</math>
|<math>4.2</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_2</math>
|<math>-11.0</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_3</math>
|<math>-41.85</math>
|<math> \frac{kJ}{mol} </math>
|-
|Solution density
|<math>\rho</math>
|<math>0.9342</math>
|<math> \frac{kg}{l} </math>
|-
|Capacity of aqueous solution
|<math>C_p</math>
|<math>3.01</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Heat transfer coefficient for cooling jacket
|<math>k_w</math>
|<math>4032</math>
|<math> \frac{kJ}{h \cdot m^2 \cdot K} </math>
|-
|Reactor surface area
|<math>A_R</math>
|<math>0.215</math>
|<math> m^2 </math>
|-
|Reactor volume
|<math>V_R</math>
|<math>10</math>
|<math> l </math>
|-
|Coolant mass
|<math>m_K</math>
|<math>5</math>
|<math> kg </math>
|-
|Capacity of coolant solution
|<math>C_{PK}</math>
|<math>2.0</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Starting concentration of subs. <math> A </math>
|<math>c_{A0}</math>
|<math>5.1</math>
|<math> \frac{mol}{l} </math>
|-
|Inflow temperature
|<math>\theta_0</math>
|<math>104.9</math>
|<math> ^\circ C </math>
|}

== Reference solution ==

"The result of a steady state optimization of the yield <math> = \frac{c_B |_S}{c_{A0}} </math> with respect to the design parameter <math> \theta_0 </math> (feed temperature) and the two controls yields the steady stae and controls"
<math> c_A = 2.1402 \frac{mol}{l}, c_B = 1.0903\frac{mol}{l}, \theta = 114.19^\circ C, \theta_K = 112.91^\circ C </math> and <math> \frac{\dot{V}}{V_R} = 14.19 h^{-1}, \dot{Q}_K = -1113.5 \frac{kJ}{h} </math>.

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Continuously Stirred Tank Reactor (TACO)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Continuously Stirred Tank Reactor problem

2019-09-02T09:51:54Z

ClemensZeile: /* Parameters */

{{Dimensions
|nd = 1
|nx = 4
|nu = 2
|nre = 2
}}The Continuously Stirred Tank Reactor problem considers a chemical reaction that produces cyclopenthenol while using up cyclepentadiene "by an acid-catalyzed electrophilic hydration in aqueous solution", an exothermal reaction which needs to be cooled. This problem can e.g. be found in <bib id="Diehl2001" />.

The inflow into the tank contains only cyclopentadiene (substance <math> A </math>) with temperature <math> \theta_0 </math> and the flow rate <math> \dot{V} </math> can be controlled. The outflow rate is the same as the inflow rate to keep the liquid level in the tank constant.
"The outflow contains a remainder of cyclopentadiene, the wanted product cyclepentenol (substance <math> B </math>) and two unwated by-products, cyclopentanediol (substance <math> C </math>) and dicyclopentadiene (substance <math> D </math>) with concentrations <math> c_A, c_B, c_C, c_D </math>." The latter two are not tracked in the problem as the substances are not of use.
The reaction scheme is given as:

<math>
\begin{array}{ccccc}
A & \overset{k_1}\rightarrow & B \overset{k_2}\rightarrow & C\\
2A & \overset{k_3}\rightarrow & D
\end{array}
</math>

where the reaction rates <math> k_i </math> are a function of the reactor temperature <math> \theta </math> via an Arrhenius law
<math> k_i(\theta) = k_{i0} \cdot \exp ( \frac{E_i}{\theta / ^\circ C + 273.15} ), \quad i=1,2,3. </math>

"The temperature <math> \theta_K </math> in the cooling jacket is held down by an external heat exchanger whose heat removal rate <math> \dot{Q}_K </math> can be controlled."

== Mathematical formulation ==
The problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{\dot{V}, \dot{Q}_K} & c_B & & \text{ at the end of reaction} \\[1.5ex]
\mbox{s.t.}
& \dot{c_A} & = & \frac{\dot{V}}{V_R} (c_{A0} - c_A) - k_1 c_A - k_3 c_A^2, \\[0.5cm]
& \dot{c_B} & = & -\frac{\dot{V}}{V_R} c_B + k_1 c_A - k_2 c_B, \\[0.5cm]
& \dot{\theta} & = & \frac{\dot{V}}{V_R} ( \theta_0 - \theta) + \frac{k_w A_R}{\rho C_p V_R} (\theta_K - \theta) - \frac{1}{\rho C_p} (k_1 c_A H_1 + k_2 c_B H_2 + k_3 c_A^2 H_3), \\[0.5cm]
& \dot{\theta_K} & = & \frac{1}{m_K C_{PK}} ( \dot{Q}_K + k_w A_R (\theta - \theta_K)),\\[0.7cm]
& c_A(0) & = & c_{A0},\\
& c_B(0) & = & 0.
\end{array}
</math>

where the various values are given in the Parameters section.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Arrhenius coefficient
|<math>k_{10}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{20}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{30}</math>
|<math> 9.043 \cdot 10^{19} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>E_1</math>
|<math> -9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_2</math>
|<math>-9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_3</math>
|<math>-8560</math>
|[-]
|-
|Reaction enthalpy
|<math>H_1</math>
|<math>4.2</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_2</math>
|<math>-11.0</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_3</math>
|<math>-41.85</math>
|<math> \frac{kJ}{mol} </math>
|-
|Solution density
|<math>\rho</math>
|<math>0.9342</math>
|<math> \frac{kg}{l} </math>
|-
|Capacity of aqueous solution
|<math>C_p</math>
|<math>3.01</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Heat transfer coefficient for cooling jacket
|<math>k_w</math>
|<math>4032</math>
|<math> \frac{kJ}{h \cdot m^2 \cdot K} </math>
|-
|Reactor surface area
|<math>A_R</math>
|<math>0.215</math>
|<math> m^2 </math>
|-
|Reactor volume
|<math>V_R</math>
|<math>10</math>
|<math> l </math>
|-
|Coolant mass
|<math>m_K</math>
|<math>5</math>
|<math> kg </math>
|-
|Capacity of coolant solution
|<math>C_{PK}</math>
|<math>2.0</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Starting concentration of subs. <math> A </math>
|<math>c_{A0}</math>
|<math>5.1</math>
|<math> \frac{mol}{l} </math>
|-
|Inflow temperature
|<math>\theta_0</math>
|<math>104.9</math>
|<math> ^\circ C </math>
|}

== Reference solution ==

"The result of a steady state optimization of the yield <math> = \frac{c_B |_S}{c_{A0}} </math> with respect to the design parameter <math> \theta_0 </math> (feed temperature) and the two controls yields the steady stae and controls"
<math> c_A = 2.1402 \frac{mol}{l}, c_B = 1.0903\frac{mol}{l}, \theta = 114.19^\circ C, \theta_K = 112.91^\circ C </math> and <math> \frac{\dot{V}}{V_R} = 14.19 h^{-1}, \dot{Q}_K = -1113.5 \frac{kJ}{h} </math>.

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Continuously Stirred Tank Reactor (TACO)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Continuously Stirred Tank Reactor problem

2019-09-02T09:51:39Z

ClemensZeile: /* Parameters */

{{Dimensions
|nd = 1
|nx = 4
|nu = 2
|nre = 2
}}The Continuously Stirred Tank Reactor problem considers a chemical reaction that produces cyclopenthenol while using up cyclepentadiene "by an acid-catalyzed electrophilic hydration in aqueous solution", an exothermal reaction which needs to be cooled. This problem can e.g. be found in <bib id="Diehl2001" />.

The inflow into the tank contains only cyclopentadiene (substance <math> A </math>) with temperature <math> \theta_0 </math> and the flow rate <math> \dot{V} </math> can be controlled. The outflow rate is the same as the inflow rate to keep the liquid level in the tank constant.
"The outflow contains a remainder of cyclopentadiene, the wanted product cyclepentenol (substance <math> B </math>) and two unwated by-products, cyclopentanediol (substance <math> C </math>) and dicyclopentadiene (substance <math> D </math>) with concentrations <math> c_A, c_B, c_C, c_D </math>." The latter two are not tracked in the problem as the substances are not of use.
The reaction scheme is given as:

<math>
\begin{array}{ccccc}
A & \overset{k_1}\rightarrow & B \overset{k_2}\rightarrow & C\\
2A & \overset{k_3}\rightarrow & D
\end{array}
</math>

where the reaction rates <math> k_i </math> are a function of the reactor temperature <math> \theta </math> via an Arrhenius law
<math> k_i(\theta) = k_{i0} \cdot \exp ( \frac{E_i}{\theta / ^\circ C + 273.15} ), \quad i=1,2,3. </math>

"The temperature <math> \theta_K </math> in the cooling jacket is held down by an external heat exchanger whose heat removal rate <math> \dot{Q}_K </math> can be controlled."

== Mathematical formulation ==
The problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{\dot{V}, \dot{Q}_K} & c_B & & \text{ at the end of reaction} \\[1.5ex]
\mbox{s.t.}
& \dot{c_A} & = & \frac{\dot{V}}{V_R} (c_{A0} - c_A) - k_1 c_A - k_3 c_A^2, \\[0.5cm]
& \dot{c_B} & = & -\frac{\dot{V}}{V_R} c_B + k_1 c_A - k_2 c_B, \\[0.5cm]
& \dot{\theta} & = & \frac{\dot{V}}{V_R} ( \theta_0 - \theta) + \frac{k_w A_R}{\rho C_p V_R} (\theta_K - \theta) - \frac{1}{\rho C_p} (k_1 c_A H_1 + k_2 c_B H_2 + k_3 c_A^2 H_3), \\[0.5cm]
& \dot{\theta_K} & = & \frac{1}{m_K C_{PK}} ( \dot{Q}_K + k_w A_R (\theta - \theta_K)),\\[0.7cm]
& c_A(0) & = & c_{A0},\\
& c_B(0) & = & 0.
\end{array}
</math>

where the various values are given in the Parameters section.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Arrhenius coefficient
|<math>k_{10}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{20}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{30}</math>
|<math> 9.043 \cdot 10^{11} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>E_1</math>
|<math> -9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_2</math>
|<math>-9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_3</math>
|<math>-8560</math>
|[-]
|-
|Reaction enthalpy
|<math>H_1</math>
|<math>4.2</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_2</math>
|<math>-11.0</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_3</math>
|<math>-41.85</math>
|<math> \frac{kJ}{mol} </math>
|-
|Solution density
|<math>\rho</math>
|<math>0.9342</math>
|<math> \frac{kg}{l} </math>
|-
|Capacity of aqueous solution
|<math>C_p</math>
|<math>3.01</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Heat transfer coefficient for cooling jacket
|<math>k_w</math>
|<math>4032</math>
|<math> \frac{kJ}{h \cdot m^2 \cdot K} </math>
|-
|Reactor surface area
|<math>A_R</math>
|<math>0.215</math>
|<math> m^2 </math>
|-
|Reactor volume
|<math>V_R</math>
|<math>10</math>
|<math> l </math>
|-
|Coolant mass
|<math>m_K</math>
|<math>5</math>
|<math> kg </math>
|-
|Capacity of coolant solution
|<math>C_{PK}</math>
|<math>2.0</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Starting concentration of subs. <math> A </math>
|<math>c_{A0}</math>
|<math>5.1</math>
|<math> \frac{mol}{l} </math>
|-
|Inflow temperature
|<math>\theta_0</math>
|<math>104.9</math>
|<math> ^\circ C </math>
|}

== Reference solution ==

"The result of a steady state optimization of the yield <math> = \frac{c_B |_S}{c_{A0}} </math> with respect to the design parameter <math> \theta_0 </math> (feed temperature) and the two controls yields the steady stae and controls"
<math> c_A = 2.1402 \frac{mol}{l}, c_B = 1.0903\frac{mol}{l}, \theta = 114.19^\circ C, \theta_K = 112.91^\circ C </math> and <math> \frac{\dot{V}}{V_R} = 14.19 h^{-1}, \dot{Q}_K = -1113.5 \frac{kJ}{h} </math>.

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Continuously Stirred Tank Reactor (TACO)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Continuously Stirred Tank Reactor problem

2019-09-02T09:51:24Z

ClemensZeile: /* Parameters */

{{Dimensions
|nd = 1
|nx = 4
|nu = 2
|nre = 2
}}The Continuously Stirred Tank Reactor problem considers a chemical reaction that produces cyclopenthenol while using up cyclepentadiene "by an acid-catalyzed electrophilic hydration in aqueous solution", an exothermal reaction which needs to be cooled. This problem can e.g. be found in <bib id="Diehl2001" />.

The inflow into the tank contains only cyclopentadiene (substance <math> A </math>) with temperature <math> \theta_0 </math> and the flow rate <math> \dot{V} </math> can be controlled. The outflow rate is the same as the inflow rate to keep the liquid level in the tank constant.
"The outflow contains a remainder of cyclopentadiene, the wanted product cyclepentenol (substance <math> B </math>) and two unwated by-products, cyclopentanediol (substance <math> C </math>) and dicyclopentadiene (substance <math> D </math>) with concentrations <math> c_A, c_B, c_C, c_D </math>." The latter two are not tracked in the problem as the substances are not of use.
The reaction scheme is given as:

<math>
\begin{array}{ccccc}
A & \overset{k_1}\rightarrow & B \overset{k_2}\rightarrow & C\\
2A & \overset{k_3}\rightarrow & D
\end{array}
</math>

where the reaction rates <math> k_i </math> are a function of the reactor temperature <math> \theta </math> via an Arrhenius law
<math> k_i(\theta) = k_{i0} \cdot \exp ( \frac{E_i}{\theta / ^\circ C + 273.15} ), \quad i=1,2,3. </math>

"The temperature <math> \theta_K </math> in the cooling jacket is held down by an external heat exchanger whose heat removal rate <math> \dot{Q}_K </math> can be controlled."

== Mathematical formulation ==
The problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{\dot{V}, \dot{Q}_K} & c_B & & \text{ at the end of reaction} \\[1.5ex]
\mbox{s.t.}
& \dot{c_A} & = & \frac{\dot{V}}{V_R} (c_{A0} - c_A) - k_1 c_A - k_3 c_A^2, \\[0.5cm]
& \dot{c_B} & = & -\frac{\dot{V}}{V_R} c_B + k_1 c_A - k_2 c_B, \\[0.5cm]
& \dot{\theta} & = & \frac{\dot{V}}{V_R} ( \theta_0 - \theta) + \frac{k_w A_R}{\rho C_p V_R} (\theta_K - \theta) - \frac{1}{\rho C_p} (k_1 c_A H_1 + k_2 c_B H_2 + k_3 c_A^2 H_3), \\[0.5cm]
& \dot{\theta_K} & = & \frac{1}{m_K C_{PK}} ( \dot{Q}_K + k_w A_R (\theta - \theta_K)),\\[0.7cm]
& c_A(0) & = & c_{A0},\\
& c_B(0) & = & 0.
\end{array}
</math>

where the various values are given in the Parameters section.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Arrhenius coefficient
|<math>k_{10}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{20}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{30}</math>
|<math> 9.043 \cdot 10^{8} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>E_1</math>
|<math> -9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_2</math>
|<math>-9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_3</math>
|<math>-8560</math>
|[-]
|-
|Reaction enthalpy
|<math>H_1</math>
|<math>4.2</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_2</math>
|<math>-11.0</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_3</math>
|<math>-41.85</math>
|<math> \frac{kJ}{mol} </math>
|-
|Solution density
|<math>\rho</math>
|<math>0.9342</math>
|<math> \frac{kg}{l} </math>
|-
|Capacity of aqueous solution
|<math>C_p</math>
|<math>3.01</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Heat transfer coefficient for cooling jacket
|<math>k_w</math>
|<math>4032</math>
|<math> \frac{kJ}{h \cdot m^2 \cdot K} </math>
|-
|Reactor surface area
|<math>A_R</math>
|<math>0.215</math>
|<math> m^2 </math>
|-
|Reactor volume
|<math>V_R</math>
|<math>10</math>
|<math> l </math>
|-
|Coolant mass
|<math>m_K</math>
|<math>5</math>
|<math> kg </math>
|-
|Capacity of coolant solution
|<math>C_{PK}</math>
|<math>2.0</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Starting concentration of subs. <math> A </math>
|<math>c_{A0}</math>
|<math>5.1</math>
|<math> \frac{mol}{l} </math>
|-
|Inflow temperature
|<math>\theta_0</math>
|<math>104.9</math>
|<math> ^\circ C </math>
|}

== Reference solution ==

"The result of a steady state optimization of the yield <math> = \frac{c_B |_S}{c_{A0}} </math> with respect to the design parameter <math> \theta_0 </math> (feed temperature) and the two controls yields the steady stae and controls"
<math> c_A = 2.1402 \frac{mol}{l}, c_B = 1.0903\frac{mol}{l}, \theta = 114.19^\circ C, \theta_K = 112.91^\circ C </math> and <math> \frac{\dot{V}}{V_R} = 14.19 h^{-1}, \dot{Q}_K = -1113.5 \frac{kJ}{h} </math>.

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Continuously Stirred Tank Reactor (TACO)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Continuously Stirred Tank Reactor problem

2019-09-02T09:51:06Z

ClemensZeile: /* Parameters */

{{Dimensions
|nd = 1
|nx = 4
|nu = 2
|nre = 2
}}The Continuously Stirred Tank Reactor problem considers a chemical reaction that produces cyclopenthenol while using up cyclepentadiene "by an acid-catalyzed electrophilic hydration in aqueous solution", an exothermal reaction which needs to be cooled. This problem can e.g. be found in <bib id="Diehl2001" />.

The inflow into the tank contains only cyclopentadiene (substance <math> A </math>) with temperature <math> \theta_0 </math> and the flow rate <math> \dot{V} </math> can be controlled. The outflow rate is the same as the inflow rate to keep the liquid level in the tank constant.
"The outflow contains a remainder of cyclopentadiene, the wanted product cyclepentenol (substance <math> B </math>) and two unwated by-products, cyclopentanediol (substance <math> C </math>) and dicyclopentadiene (substance <math> D </math>) with concentrations <math> c_A, c_B, c_C, c_D </math>." The latter two are not tracked in the problem as the substances are not of use.
The reaction scheme is given as:

<math>
\begin{array}{ccccc}
A & \overset{k_1}\rightarrow & B \overset{k_2}\rightarrow & C\\
2A & \overset{k_3}\rightarrow & D
\end{array}
</math>

where the reaction rates <math> k_i </math> are a function of the reactor temperature <math> \theta </math> via an Arrhenius law
<math> k_i(\theta) = k_{i0} \cdot \exp ( \frac{E_i}{\theta / ^\circ C + 273.15} ), \quad i=1,2,3. </math>

"The temperature <math> \theta_K </math> in the cooling jacket is held down by an external heat exchanger whose heat removal rate <math> \dot{Q}_K </math> can be controlled."

== Mathematical formulation ==
The problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{\dot{V}, \dot{Q}_K} & c_B & & \text{ at the end of reaction} \\[1.5ex]
\mbox{s.t.}
& \dot{c_A} & = & \frac{\dot{V}}{V_R} (c_{A0} - c_A) - k_1 c_A - k_3 c_A^2, \\[0.5cm]
& \dot{c_B} & = & -\frac{\dot{V}}{V_R} c_B + k_1 c_A - k_2 c_B, \\[0.5cm]
& \dot{\theta} & = & \frac{\dot{V}}{V_R} ( \theta_0 - \theta) + \frac{k_w A_R}{\rho C_p V_R} (\theta_K - \theta) - \frac{1}{\rho C_p} (k_1 c_A H_1 + k_2 c_B H_2 + k_3 c_A^2 H_3), \\[0.5cm]
& \dot{\theta_K} & = & \frac{1}{m_K C_{PK}} ( \dot{Q}_K + k_w A_R (\theta - \theta_K)),\\[0.7cm]
& c_A(0) & = & c_{A0},\\
& c_B(0) & = & 0.
\end{array}
</math>

where the various values are given in the Parameters section.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Arrhenius coefficient
|<math>k_{10}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{20}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{30}</math>
|<math> 9.043 \cdot 10^8 </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>E_1</math>
|<math> -9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_2</math>
|<math>-9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_3</math>
|<math>-8560</math>
|[-]
|-
|Reaction enthalpy
|<math>H_1</math>
|<math>4.2</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_2</math>
|<math>-11.0</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_3</math>
|<math>-41.85</math>
|<math> \frac{kJ}{mol} </math>
|-
|Solution density
|<math>\rho</math>
|<math>0.9342</math>
|<math> \frac{kg}{l} </math>
|-
|Capacity of aqueous solution
|<math>C_p</math>
|<math>3.01</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Heat transfer coefficient for cooling jacket
|<math>k_w</math>
|<math>4032</math>
|<math> \frac{kJ}{h \cdot m^2 \cdot K} </math>
|-
|Reactor surface area
|<math>A_R</math>
|<math>0.215</math>
|<math> m^2 </math>
|-
|Reactor volume
|<math>V_R</math>
|<math>10</math>
|<math> l </math>
|-
|Coolant mass
|<math>m_K</math>
|<math>5</math>
|<math> kg </math>
|-
|Capacity of coolant solution
|<math>C_{PK}</math>
|<math>2.0</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Starting concentration of subs. <math> A </math>
|<math>c_{A0}</math>
|<math>5.1</math>
|<math> \frac{mol}{l} </math>
|-
|Inflow temperature
|<math>\theta_0</math>
|<math>104.9</math>
|<math> ^\circ C </math>
|}

== Reference solution ==

"The result of a steady state optimization of the yield <math> = \frac{c_B |_S}{c_{A0}} </math> with respect to the design parameter <math> \theta_0 </math> (feed temperature) and the two controls yields the steady stae and controls"
<math> c_A = 2.1402 \frac{mol}{l}, c_B = 1.0903\frac{mol}{l}, \theta = 114.19^\circ C, \theta_K = 112.91^\circ C </math> and <math> \frac{\dot{V}}{V_R} = 14.19 h^{-1}, \dot{Q}_K = -1113.5 \frac{kJ}{h} </math>.

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Continuously Stirred Tank Reactor (TACO)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Continuously Stirred Tank Reactor problem

2019-09-02T09:50:11Z

ClemensZeile: /* Parameters */

{{Dimensions
|nd = 1
|nx = 4
|nu = 2
|nre = 2
}}The Continuously Stirred Tank Reactor problem considers a chemical reaction that produces cyclopenthenol while using up cyclepentadiene "by an acid-catalyzed electrophilic hydration in aqueous solution", an exothermal reaction which needs to be cooled. This problem can e.g. be found in <bib id="Diehl2001" />.

The inflow into the tank contains only cyclopentadiene (substance <math> A </math>) with temperature <math> \theta_0 </math> and the flow rate <math> \dot{V} </math> can be controlled. The outflow rate is the same as the inflow rate to keep the liquid level in the tank constant.
"The outflow contains a remainder of cyclopentadiene, the wanted product cyclepentenol (substance <math> B </math>) and two unwated by-products, cyclopentanediol (substance <math> C </math>) and dicyclopentadiene (substance <math> D </math>) with concentrations <math> c_A, c_B, c_C, c_D </math>." The latter two are not tracked in the problem as the substances are not of use.
The reaction scheme is given as:

<math>
\begin{array}{ccccc}
A & \overset{k_1}\rightarrow & B \overset{k_2}\rightarrow & C\\
2A & \overset{k_3}\rightarrow & D
\end{array}
</math>

where the reaction rates <math> k_i </math> are a function of the reactor temperature <math> \theta </math> via an Arrhenius law
<math> k_i(\theta) = k_{i0} \cdot \exp ( \frac{E_i}{\theta / ^\circ C + 273.15} ), \quad i=1,2,3. </math>

"The temperature <math> \theta_K </math> in the cooling jacket is held down by an external heat exchanger whose heat removal rate <math> \dot{Q}_K </math> can be controlled."

== Mathematical formulation ==
The problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{\dot{V}, \dot{Q}_K} & c_B & & \text{ at the end of reaction} \\[1.5ex]
\mbox{s.t.}
& \dot{c_A} & = & \frac{\dot{V}}{V_R} (c_{A0} - c_A) - k_1 c_A - k_3 c_A^2, \\[0.5cm]
& \dot{c_B} & = & -\frac{\dot{V}}{V_R} c_B + k_1 c_A - k_2 c_B, \\[0.5cm]
& \dot{\theta} & = & \frac{\dot{V}}{V_R} ( \theta_0 - \theta) + \frac{k_w A_R}{\rho C_p V_R} (\theta_K - \theta) - \frac{1}{\rho C_p} (k_1 c_A H_1 + k_2 c_B H_2 + k_3 c_A^2 H_3), \\[0.5cm]
& \dot{\theta_K} & = & \frac{1}{m_K C_{PK}} ( \dot{Q}_K + k_w A_R (\theta - \theta_K)),\\[0.7cm]
& c_A(0) & = & c_{A0},\\
& c_B(0) & = & 0.
\end{array}
</math>

where the various values are given in the Parameters section.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Arrhenius coefficient
|<math>k_{10}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{20}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{30}</math>
|<math> 9.043 \cdot 10^9 </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>E_1</math>
|<math> -9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_2</math>
|<math>-9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_3</math>
|<math>-8560</math>
|[-]
|-
|Reaction enthalpy
|<math>H_1</math>
|<math>4.2</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_2</math>
|<math>-11.0</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_3</math>
|<math>-41.85</math>
|<math> \frac{kJ}{mol} </math>
|-
|Solution density
|<math>\rho</math>
|<math>0.9342</math>
|<math> \frac{kg}{l} </math>
|-
|Capacity of aqueous solution
|<math>C_p</math>
|<math>3.01</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Heat transfer coefficient for cooling jacket
|<math>k_w</math>
|<math>4032</math>
|<math> \frac{kJ}{h \cdot m^2 \cdot K} </math>
|-
|Reactor surface area
|<math>A_R</math>
|<math>0.215</math>
|<math> m^2 </math>
|-
|Reactor volume
|<math>V_R</math>
|<math>10</math>
|<math> l </math>
|-
|Coolant mass
|<math>m_K</math>
|<math>5</math>
|<math> kg </math>
|-
|Capacity of coolant solution
|<math>C_{PK}</math>
|<math>2.0</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Starting concentration of subs. <math> A </math>
|<math>c_{A0}</math>
|<math>5.1</math>
|<math> \frac{mol}{l} </math>
|-
|Inflow temperature
|<math>\theta_0</math>
|<math>104.9</math>
|<math> ^\circ C </math>
|}

== Reference solution ==

"The result of a steady state optimization of the yield <math> = \frac{c_B |_S}{c_{A0}} </math> with respect to the design parameter <math> \theta_0 </math> (feed temperature) and the two controls yields the steady stae and controls"
<math> c_A = 2.1402 \frac{mol}{l}, c_B = 1.0903\frac{mol}{l}, \theta = 114.19^\circ C, \theta_K = 112.91^\circ C </math> and <math> \frac{\dot{V}}{V_R} = 14.19 h^{-1}, \dot{Q}_K = -1113.5 \frac{kJ}{h} </math>.

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Continuously Stirred Tank Reactor (TACO)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Continuously Stirred Tank Reactor problem

2019-09-02T09:49:52Z

ClemensZeile: /* Parameters */

Continuously Stirred Tank Reactor problem

2019-09-02T09:49:34Z

ClemensZeile: /* Parameters */

{{Dimensions
|nd = 1
|nx = 4
|nu = 2
|nre = 2
}}The Continuously Stirred Tank Reactor problem considers a chemical reaction that produces cyclopenthenol while using up cyclepentadiene "by an acid-catalyzed electrophilic hydration in aqueous solution", an exothermal reaction which needs to be cooled. This problem can e.g. be found in <bib id="Diehl2001" />.

The inflow into the tank contains only cyclopentadiene (substance <math> A </math>) with temperature <math> \theta_0 </math> and the flow rate <math> \dot{V} </math> can be controlled. The outflow rate is the same as the inflow rate to keep the liquid level in the tank constant.
"The outflow contains a remainder of cyclopentadiene, the wanted product cyclepentenol (substance <math> B </math>) and two unwated by-products, cyclopentanediol (substance <math> C </math>) and dicyclopentadiene (substance <math> D </math>) with concentrations <math> c_A, c_B, c_C, c_D </math>." The latter two are not tracked in the problem as the substances are not of use.
The reaction scheme is given as:

<math>
\begin{array}{ccccc}
A & \overset{k_1}\rightarrow & B \overset{k_2}\rightarrow & C\\
2A & \overset{k_3}\rightarrow & D
\end{array}
</math>

where the reaction rates <math> k_i </math> are a function of the reactor temperature <math> \theta </math> via an Arrhenius law
<math> k_i(\theta) = k_{i0} \cdot \exp ( \frac{E_i}{\theta / ^\circ C + 273.15} ), \quad i=1,2,3. </math>

"The temperature <math> \theta_K </math> in the cooling jacket is held down by an external heat exchanger whose heat removal rate <math> \dot{Q}_K </math> can be controlled."

== Mathematical formulation ==
The problem is given by


<math>
\begin{array}{llcl}
\displaystyle \max_{\dot{V}, \dot{Q}_K} & c_B & & \text{ at the end of reaction} \\[1.5ex]
\mbox{s.t.}
& \dot{c_A} & = & \frac{\dot{V}}{V_R} (c_{A0} - c_A) - k_1 c_A - k_3 c_A^2, \\[0.5cm]
& \dot{c_B} & = & -\frac{\dot{V}}{V_R} c_B + k_1 c_A - k_2 c_B, \\[0.5cm]
& \dot{\theta} & = & \frac{\dot{V}}{V_R} ( \theta_0 - \theta) + \frac{k_w A_R}{\rho C_p V_R} (\theta_K - \theta) - \frac{1}{\rho C_p} (k_1 c_A H_1 + k_2 c_B H_2 + k_3 c_A^2 H_3), \\[0.5cm]
& \dot{\theta_K} & = & \frac{1}{m_K C_{PK}} ( \dot{Q}_K + k_w A_R (\theta - \theta_K)),\\[0.7cm]
& c_A(0) & = & c_{A0},\\
& c_B(0) & = & 0.
\end{array}
</math>

where the various values are given in the Parameters section.

== Parameters ==
These fixed values are used within the model.
{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Arrhenius coefficient
|<math>k_{10}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{20}</math>
|<math> 1.287 \cdot 10^{12} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>k_{30}</math>
|<math> 9.043 \cdot 10^{19} </math>
|<math> h^{-1} </math>
|-
|Arrhenius coefficient
|<math>E_1</math>
|<math> -9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_2</math>
|<math>-9758.3</math>
|[-]
|-
|Arrhenius coefficient
|<math>E_3</math>
|<math>-8560</math>
|[-]
|-
|Reaction enthalpy
|<math>H_1</math>
|<math>4.2</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_2</math>
|<math>-11.0</math>
|<math> \frac{kJ}{mol} </math>
|-
|Reaction enthalpy
|<math>H_3</math>
|<math>-41.85</math>
|<math> \frac{kJ}{mol} </math>
|-
|Solution density
|<math>\rho</math>
|<math>0.9342</math>
|<math> \frac{kg}{l} </math>
|-
|Capacity of aqueous solution
|<math>C_p</math>
|<math>3.01</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Heat transfer coefficient for cooling jacket
|<math>k_w</math>
|<math>4032</math>
|<math> \frac{kJ}{h \cdot m^2 \cdot K} </math>
|-
|Reactor surface area
|<math>A_R</math>
|<math>0.215</math>
|<math> m^2 </math>
|-
|Reactor volume
|<math>V_R</math>
|<math>10</math>
|<math> l </math>
|-
|Coolant mass
|<math>m_K</math>
|<math>5</math>
|<math> kg </math>
|-
|Capacity of coolant solution
|<math>C_{PK}</math>
|<math>2.0</math>
|<math> \frac{kJ}{kg \cdot K} </math>
|-
|Starting concentration of subs. <math> A </math>
|<math>c_{A0}</math>
|<math>5.1</math>
|<math> \frac{mol}{l} </math>
|-
|Inflow temperature
|<math>\theta_0</math>
|<math>104.9</math>
|<math> ^\circ C </math>
|}

== Reference solution ==

"The result of a steady state optimization of the yield <math> = \frac{c_B |_S}{c_{A0}} </math> with respect to the design parameter <math> \theta_0 </math> (feed temperature) and the two controls yields the steady stae and controls"
<math> c_A = 2.1402 \frac{mol}{l}, c_B = 1.0903\frac{mol}{l}, \theta = 114.19^\circ C, \theta_K = 112.91^\circ C </math> and <math> \frac{\dot{V}}{V_R} = 14.19 h^{-1}, \dot{Q}_K = -1113.5 \frac{kJ}{h} </math>.

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Continuously Stirred Tank Reactor (TACO)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

Catalyst mixing problem

2018-01-12T19:15:36Z

ClemensZeile:

{{Dimensions
|nd = 1
|nx = 2
|nu = 1
|nc = 2
|nre = 2
}}The Catalyst mixing problem seeks an optimal policy for mixing two catalysts "along the length of a tubular plug
ow reactor involving several reactions". (Cite and problem taken from the [http://www.mcs.anl.gov/~more/cops/ COPS library])

== Mathematical formulation ==

The problem is given by


<math>
\begin{array}{llcl}
\displaystyle \min_{x, w} &-1 + x_1(t_f) + x_2(t_f) \\[1.5ex]
\mbox{s.t.}
& \dot{x}_1 & = & w(t) ( 10 x_2(t) - x_1(t)), \\
& \dot{x}_2 & = & w(t) ( x_1(t) - 10 x_2(t)) - (1 - w(t)) \, x_2(t) , \\
& x(t_0) &=& (1, 0)^T, \\
& w(t) &\in& \{0,1\}.
\end{array}
</math>


== Parameters ==

In this model the parameters used are <math> t_0 = 0, \, \, t_f = 1 </math>.

== Reference Solution ==
If the problem is relaxed, i.e., we demand that w(t) be in the continuous interval [0, 1] instead of the binary choice \{0,1\}, the optimal solution can be determined by means of direct optimal control.

<gallery caption="Reference solution plots" widths="180px" heights="140px" perrow="2">
Image:Catalyst_Mixing_Problem_Performance.png| Results with relaxed controls and collocation from the [http://www.mcs.anl.gov/~more/cops/ COPS library]
Image:Catalyst Mixing Controls.png| Optimal relaxed controls showing a bang-bang structure.
</gallery>

== Source Code ==

Model descriptions are available in

* [[:Category:AMPL/TACO | AMPL/TACO code]] at [[Catalyst mixing problem (TACO)]]


[[Category:MIOCP]]
[[Category:ODE model]]
[[Category:Chemical engineering]]

File:Catalyst Mixing Controls.png

2018-01-12T19:14:49Z

ClemensZeile:

File:Catalyst Mixing Problem Performance.png

2018-01-12T19:13:01Z

ClemensZeile: