Difference between revisions of "Bioreactor (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 = 0*ones(N,1);
 input.yN = 0;
 
 input.W = D;
 input.WN = D;
 
 %% 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')