Bioreactor (ACADO)
From mintOC
Revision as of 20:23, 26 December 2015 by FelixMueller (Talk | contribs) (Created page with "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 m...")
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 = 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')