Lotka Volterra fishing problem (ACADO)

This pages states the code to solve the (ralaxed) Lotka Volterra fishing problem with the code generation fo the ACADO Toolkit from Matlab (via Matlab interface). Hence one needs Matlab and ACADO. The m-file should be located in the matlab interface of ACADO and the path <ACADOROOT>/interfaces/matlab/examples/codegeneratio/nmpc/. While we use the NMPC more precisely the RTI capability of the ACADO Code generation, this features are misused to implement a 20 step SQP solver without stopping criteria.

ACADO

The source code for a fixed grid discretisation with fixed stepsize integrator.

% implements the lotka volterra fishing problem on mintoc.de
 
clc;
clear all;
close all;
 
Ts = 0.1;
EXPORT = 1;
 
%% Variables
DifferentialState x0 x1;
Control w;
 
n_XD = length(diffStates);
n_U = length(controls);
 
% Constants
c0 = 0.4;
c1 = 0.2;
 
%% Differential Equation
 
f = dot([x0; x1]) == [ x0-x0*x1-c0*x0*w; ...
    -x1 + x0*x1-c1*x1*w];
 
h = [diffStates];
hN = [diffStates];
 
 
 
%% MPCexport
acadoSet('problemname', 'mpc');
 
N = 200;
ocp = acado.OCP( 0.0, 12, N );
 
W_mat = eye(n_XD,n_XD);
WN_mat = eye(n_XD,n_XD);
W = acado.BMatrix(W_mat);
WN = acado.BMatrix(WN_mat);
 
ocp.minimizeLSQ( W, h );
ocp.minimizeLSQEndTerm( WN, hN );
 
ocp.subjectTo( 0 <= w <= 1 );
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',        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
 
%% PARAMETERS OPTIMIZATION
X0 = [0.5 0.7];
input.x0=X0';
Xref = [1 1];
input.x = repmat(Xref,N+1,1);
Xref = repmat(Xref,N,1);
input.od = [];
 
Uref = zeros(N,n_U);
input.u = Uref;
 
input.y = [Xref(1:N,:)];
input.yN = Xref(N,:).';
 
input.W = diag([1 1]);
input.WN = diag([1 1]);
 
%% SOLVER LOOP
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
 
t_end = 12;
 
% States
figure(1)
plot([0:t_end/N:t_end],output.x)
ylabel('States')
xlabel('Time')
legend('x1','x2')
 
% Control
figure(2)
plot([0:t_end/N:t_end-t_end/N],output.u)
ylabel('Control (u)')
xlabel('Time')