The ACADO manual

Houska, Boris, et al. "ACADO toolkit user’s manual." (2009): 3.

describes the toolkit as follows:

"
ACADO Toolkit is a software environment and algorithm collection written in C++ for automatic control and dynamic optimization. It provides a general framework for using a great variety of algorithms for direct optimal control, including model predictive control as well as state and parameter estimation. It also provides (stand-alone) efficiently implemented Runge-Kutta and BDF integrators for the simulation of ODE’s and DAE’s.
"

There a list of key properties is also given. The toolkit is open source, should be user-friendly. The code should be easy to extend (i.e. linking to it is simple) and self contained (no external software is required). Due to this properties the ACADO toolkit has become quite popular in research with Optimal Control Problems.

[[Category: Implementation]]

Category:ACADO

2016-01-28T10:12:22Z

TobiasWeber:

2015-12-09T14:53:09Z

TobiasWeber: /* Source Code */

The bioreactor example is an easy bioreactor with an substrate that is converted to a product by the biomass in the reactor. It has three states and a control that is describing the feed concentration of substrate. It is taken from the examples folder of the ACADO toolkit described in:

Houska, Boris, Hans Joachim Ferreau, and Moritz Diehl. "ACADO toolkit—An open‐source framework for automatic control and dynamic optimization."
Optimal Control Applications and Methods 32.3 (2011): 298-312.

Originally the problem seems to be motivated by:

VERSYCK, KARINA J., and JAN F. VAN IMPE. "Feed rate optimization for fed-batch bioreactors: From optimal process performance to optimal parameter estimation."
Chemical Engineering Communications 172.1 (1999): 107-124.

== Model Formulation ==

The dynamic model is an [[:Category:ODE model|ODE model]]:


<math>
\begin{array}{rcl}
\dot{X}&=&-DX+\mu X \\
\dot{S}&=& D(S_{f}-S)-\mu /Y_{xs} X \\
\dot{P}&=&-DP+ (\alpha \mu +\beta) X.
\end{array}
</math>


The three states describe the concentration of the biomass (<math>X</math>), the substrate (<math>S</math>), and the product (<math>P</math>) in the reactor. In steady state the feed and outlet are equal and dilute all three concentrations with a ratio <math>D</math>. The biomass grows with a rate
<math>\mu</math>, while it eats up the substrate with the rate <math>\mu/Y_{xs}</math> and produces product at a rate <math>(\alpha \mu +\beta)</math>. The rate <math>\mu</math> is given by:

<math>\mu = \mu_{m}*(1-P/P_{m})*S/(K_m+S+S^2/K_i)</math>

The fixed parameters (constants) of the model are as follows.

{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Dilution
|<math>D</math>
|0.15
|[-]
|-
|Rate coefficient
|<math>K_i</math>
|22
|[-]
|-
|Rate coefficient
|<math>K_m</math>
|1.2
|[-]
|-
|Rate coefficient
|<math>P_m</math>
|50
|[-]
|-
|Substrate to Biomass rate
|<math>Y_{xs}</math>
|0.4
|[-]
|-
|Linear slope
|<math>\alpha</math>
|2.2
|[-]
|-
|Linear intercept
|<math>\beta</math>
|0.2
|[-]
|-
|Maximal growth rate
|<math>\mu_m</math>
|0.48
|[-]
|-
|}

== Optimal Control Problem ==

Writing shortly for the states in vector notation <math>x=(X,S,P)^T</math> the OCP reads:


<math>
\begin{array}{cl}
\displaystyle \min_{x,S_f} & J(x,S_f)\\[1.5ex]
\mbox{s.t.} & \dot{x} = f(x,S_f), \forall \, t \in [0,48]\\
& x(0) = (6.5,12,22)^T \\
& x \in \R^3,\,S_f \in [28.7,40].
\end{array}
</math>


=== Objective ===

<math>
J(x,S_f)=\int_0^{48}D(S_f-P)^2dt
</math>


== Reference Solution ==

Here we present the reference solution of the reimplemented example in the ACADO code generation with matlab. The source code is given in the next section.

<gallery caption="Reference solution" widths="551px" heights="390px" perrow="1">
Image:ACADO_bioreactor.png| Optimal solution for the ACADO example.
</gallery>

== Source Code ==

This is the implementation of the ACADO bioreactor example with the matlab interface and the code generation of ACADO.

% implements the bioreactor example of ACADO for the matlab interface and
% the code generation

clc;
clear all;
close all;

Ts = 0.1;
EXPORT = 1;

%% Variables
DifferentialState X S P;
Control Sf;

n_XD = length(diffStates);
n_U = length(controls);

%% Constants
D = 0.15;
Ki = 22.0;
Km = 1.2 ;
Pm = 50.0;
Yxs = 0.4 ;
alpha = 2.2 ;
beta = 0.2 ;
mum = 0.48;
Sfmin = 28.7;
Sfmax = 40.0;
t_start = 0.0;
t_end = 48.0;
N = 20;

%% Differential Equation
mu = mum*(1-P/Pm)*S/(Km+S+S^2/Ki);

f = dot([X;S;P]) == [-D*X+mu*X;...
D*(Sf-S)-(mu/Yxs)*X;...
-D*P+(alpha*mu+beta)*X];

% output
h = P-Sf;
hN = P;

%% MPCexport
acadoSet('problemname', 'mpc');

ocp = acado.OCP( t_start, t_end, N );

W_mat = D;
WN_mat = D;
W = acado.BMatrix(W_mat);
WN = acado.BMatrix(WN_mat);

ocp.minimizeLSQ( W, h );
ocp.minimizeLSQEndTerm( WN, hN );

ocp.subjectTo( Sfmin <= Sf <= Sfmax );
ocp.setModel(f);

mpc = acado.OCPexport( ocp );
mpc.set( 'HESSIAN_APPROXIMATION', 'GAUSS_NEWTON' );
mpc.set( 'DISCRETIZATION_TYPE', 'MULTIPLE_SHOOTING' );
mpc.set( 'SPARSE_QP_SOLUTION', 'FULL_CONDENSING_N2');
mpc.set( 'INTEGRATOR_TYPE', 'INT_IRK_GL4' );
mpc.set( 'NUM_INTEGRATOR_STEPS', 10*N );
mpc.set( 'QP_SOLVER', 'QP_QPOASES' );
mpc.set( 'HOTSTART_QP', 'NO' );
mpc.set( 'LEVENBERG_MARQUARDT', 1e-10 );

if EXPORT
mpc.exportCode( 'export_MPC' );
copyfile('../../../../../../external_packages/qpoases', 'export_MPC/qpoases')

cd export_MPC
make_acado_solver('../acado_MPCstep')
cd ..
end

%% CONSTANTS FOR OPTIMIZATION
X0 = [6.5 12.0 22.0];
input.x0=X0';
Xref = [0 0 0];
input.x = repmat(Xref,N+1,1);
Xref = repmat(Xref,N,1);
input.od = [];

Uref = zeros(N,n_U);
input.u = Uref;

input.y = Sfmin*ones(N,1);
input.yN = Sfmin;

input.W = 1;
input.WN = 1;

%% SOLVER LOOP (SQP - Gauss newton)
display('------------------------------------------------------------------')
display(' SOLVER Loop' )
display('------------------------------------------------------------------')

for i=1:20
tic
% Solve NMPC OCP
output = acado_MPCstep(input);

input.x=output.x;
input.u=output.u;

disp([' (RTI step: ' num2str(output.info.cpuTime*1e6) ' µs)'])

end

%% PLOT RESULTS

% States
figure(1)
plot([t_start:t_end/N:t_end],output.x)
ylabel('States')
xlabel('Time')
legend('X','S','P')

% Control
figure(2)
plot([t_start:t_end/N:t_end-t_end/N],output.u)
ylabel('Control (Sf)')
xlabel('Time')

% one figure for all
figure(3)
subplot(2,2,1)
plot([t_start:t_end/N:t_end],output.x(:,1))
ylabel('X')
xlabel('Time')
subplot(2,2,2)
plot([t_start:t_end/N:t_end],output.x(:,2))
ylabel('S')
xlabel('Time')
subplot(2,2,3)
plot([t_start:t_end/N:t_end],output.x(:,3))
ylabel('P')
xlabel('Time')
subplot(2,2,4)
plot([t_start:t_end/N:t_end-t_end/N],output.u)
ylabel('Sf')
xlabel('Time')

Bioreactor

2015-12-09T14:47:11Z

TobiasWeber: /* Source Code */

Bioreactor

2015-12-09T14:45:39Z

TobiasWeber: /* Reference Solution */

Bioreactor

2015-12-09T14:44:43Z

TobiasWeber: /* Reference Solution */

File:ACADO bioreactor.png

The bioreactor example is an easy bioreactor with an substrate that is converted to a product by the biomass in the reactor. It has three states and a control that is describing the feed concentration of substrate. It is taken from the examples folder of the ACADO toolkit described in:

Houska, Boris, Hans Joachim Ferreau, and Moritz Diehl. "ACADO toolkit—An open‐source framework for automatic control and dynamic optimization."
Optimal Control Applications and Methods 32.3 (2011): 298-312.

Originally the problem seems to be motivated by:

VERSYCK, KARINA J., and JAN F. VAN IMPE. "Feed rate optimization for fed-batch bioreactors: From optimal process performance to optimal parameter estimation."
Chemical Engineering Communications 172.1 (1999): 107-124.

== Model Formulation ==

The dynamic model is an [[:Category:ODE model|ODE model]]:


<math>
\begin{array}{rcl}
\dot{X}&=&-DX+\mu X \\
\dot{S}&=& D(S_{f}-S)-\mu /Y_{xs} X \\
\dot{P}&=&-DP+ (\alpha \mu +\beta) X.
\end{array}
</math>


The three states describe the concentration of the biomass (<math>X</math>), the substrate (<math>S</math>), and the product (<math>P</math>) in the reactor. In steady state the feed and outlet are equal and dilute all three concentrations with a ratio <math>D</math>. The biomass grows with a rate
<math>\mu</math>, while it eats up the substrate with the rate <math>\mu/Y_{xs}</math> and produces product at a rate <math>(\alpha \mu +\beta)</math>. The rate <math>\mu</math> is given by:

<math>\mu = \mu_{m}*(1-P/P_{m})*S/(K_m+S+S^2/K_i)</math>

The fixed parameters (constants) of the model are as follows.

{| class="wikitable"
|+Parameters
|-
|Name
|Symbol
|Value
|Unit
|-
|Dilution
|<math>D</math>
|0.15
|[-]
|-
|Rate coefficient
|<math>K_i</math>
|22
|[-]
|-
|Rate coefficient
|<math>K_m</math>
|1.2
|[-]
|-
|Rate coefficient
|<math>P_m</math>
|50
|[-]
|-
|Substrate to Biomass rate
|<math>Y_{xs}</math>
|0.4
|[-]
|-
|Linear slope
|<math>\alpha</math>
|2.2
|[-]
|-
|Linear intercept
|<math>\beta</math>
|0.2
|[-]
|-
|Maximal growth rate
|<math>\mu_m</math>
|0.48
|[-]
|-
|}

== Optimal Control Problem ==

In this section the Optimal Control Problem (OCP) with the dynamic model as equality constraint is specified. In general every OCP has either a Lagrange or Mayer type objective or both.
First the general structure of the OCP is given to show the type of objective (linear, quadratic, ...) and constraints (equality and/or inequality, linear or nonlinear) as well as the variables that are free to optimize.
All of the above dynamic models can be used, in the following we just write for all:


<math>
\dot{x}=f(x,u,p,t) \forall \, t \in I.
</math>


Implicit systems would be treated similarly. Then the general OCP is writte as:


<math>
\begin{array}{cl}
\displaystyle \min_{x, u, p} & J(x,u,t)\\[1.5ex]
\mbox{s.t.} & \dot{x} = f(x,u,p,t), \forall \, t \in I\\
& 0 = g(x(t_o),x(t_f),p) \\
& 0 \ge c(x,u,p), \forall \, t \in I\\
& 0 = h(x,u,p), \forall \, t \in I \\
& x \in \mathcal{X},\,u \in \mathcal{U},\, p \in P.
\end{array}
</math>


This is the most general formulation, where the initial value, parameters, and inputs are free. The states and inputs are vectors (or matrices) of functions while the parameters are just a vector of variables. If the initial (or final) state is fixed, one just adds a constraint <math>x(t_0)=x_0</math> and specifies the value(s) of <math>x_0</math> in a list or table.
If the parameters are fixed, then one just removes them from the list of optimization variables and states their values in a list (same for point conditions inside of the interval <math>I</math>). If the OCP is a parameter estimation problem without inputs,
the inputs just dont appear in the problem above at all, neither as optimization variable nor as part of the constraints. In contrast if they are just fixed this is done by the equality constraints. In this setup the above problem is very general and by
stating the different objective functions and constraint functions one can cover a wide range of OCPs.

Bioreactor

2015-12-08T14:48:13Z

TobiasWeber:

Bioreactor

2015-12-08T14:46:19Z

TobiasWeber:

Bioreactor

2015-12-08T14:42:44Z

TobiasWeber:

Bioreactor

2015-12-08T14:36:54Z

TobiasWeber: