Optimization and Control of a Fed-Batch Reactor Using Nonlinear MPC

Fed-Batch Chemical Reactor

The following irreversible and exothermic reactions occur in the batch reactor [1]:

A + B => C (desired product) C => D (undesired product)

The batch begins with the reactor partially filled with known concentrations of reactantsAandB. The batch reacts for 0.5 hours, during which additionalBcan be added and the reactor temperature can be changed.

The nonlinear model of the batch reactor is defined in thefedbatch_StateFcnandfedbatch_OutputFcnfunctions. This system has the following inputs, states, and outputs.

Manipulated Variables

Measured disturbance

States

Outputs

The goal is to maximize the production ofC(y1) at the end of the batch process. During the batch process, the following operating constraints must be satisfied:

Specify the nominal operating condition at the beginning of the batch process.

c_A0 = 10; c_B0 = 1.167; c_C0 = 0; V0 = 1; T0 = 50 + 273.15; c_Bin = 20;

Specify the nominal states.

x0 = zeros(3,1); x0(1) = c_A0*V0; x0(2) = x0(1) + c_C0*V0; x0(3) = V0; x0(4) = T0;

Specify the nominal inputs.

u0 = zeros(3,1); u0(2) = 40; u0(3) = c_Bin;

Specify the nominal outputs.

y0 = fedbatch_OutputFcn(x0,u0);

Nonlinear MPC Design to Optimize Batch Operation

Create a nonlinear MPC object with 4 states, 3 outputs, 2 manipulated variables, and 1 measured disturbance.

nlmpcobj_Plan = nlmpc(4, 3,'MV', [1,2],'MD', 3);

Zero weights are applied to one or more OVs because there are fewer MVs than OVs.

Given the expected batch durationTf, choose the controller sample timeTsand prediction horizon.

Tf = 0.5; N = 50; Ts = Tf/N; nlmpcobj_Plan.Ts = Ts; nlmpcobj_Plan.PredictionHorizon = N;

If you set the control horizon equal to the prediction horizon, there will be 50 free control moves, which leads to a total of 100 decision variables because the plant has two manipulated variables. To reduce the number of decision variables, you can specify control horizon using blocking moves. Divide the prediction horizon into 8 blocks, which represents 8 free control moves. Each of the first seven blocks lasts seven prediction steps. Doing so reduces the number of decision variables to 16.

nlmpcobj_Plan.ControlHorizon = [7 7 7 7 7 7 7 1];

指定的非线性model in the controller. The functionfedbatch_StateFcnDTconverts the continuous-time model to discrete time using a multi-step Forward Euler integration formula.

nlmpcobj_Plan.Model.StateFcn = @(x,u) fedbatch_StateFcnDT(x,u,Ts); nlmpcobj_Plan.Model.OutputFcn = @(x,u) fedbatch_OutputFcn(x,u); nlmpcobj_Plan.Model.IsContinuousTime = false;

Specify the bounds for feed rate ofB.

nlmpcobj_Plan.MV(1).Min = 0; nlmpcobj_Plan.MV(1).Max = 1;

Specify the bounds for the reactor temperature setpoint.

nlmpcobj_Plan.MV(2).Min = 20; nlmpcobj_Plan.MV(2).Max = 50;

Specify the upper bound for the heat removal rate. The true constraint is1.5e5. Since nonlinear MPC can only enforce constraints at the sampling instants, use a safety margin of0.05e5to prevent a constraint violation between sampling instants.

nlmpcobj_Plan.OV(2).Max = 1.45e5;

Specify the upper bound for the liquid volume in the reactor.

nlmpcobj_Plan.OV(3).Max = 1.1;

Since the goal is to maximizey1, the amount ofCin the reactor at the end of the batch time, specify a custom cost function that replaces the default quadratic cost. Sincey1 = x2-x1, define the custom cost to be minimized asx1-x2using an anonymous function.

nlmpcobj_Plan.Optimization.CustomCostFcn = @(X,U,e,data) X(end,1)-X(end,2); nlmpcobj_Plan.Optimization.ReplaceStandardCost = true;

To configure the manipulated variables to vary linearly with time within each block, select piecewise linear interpolation. By default, nonlinear MPC keeps manipulated variables constant within each block, using piecewise constant interpolation, which might be too restrictive for an optimal trajectory planning problem.

nlmpcobj_Plan.Optimization.MVInterpolationOrder = 1;

Use the default nonlinear programming solverfminconto solve the nonlinear MPC problem. For this example, set the solver step tolerance to help achieve first order optimality.

nlmpcobj_Plan.Optimization.SolverOptions.StepTolerance = 1e-8;

Before carrying out optimization, check whether all the custom functions satisfy NLMPC requirements using thevalidateFcnscommand.

validateFcns(nlmpcobj_Plan, x0, u0(1:2), u0(3));

Model.StateFcn is OK. Model.OutputFcn is OK. Optimization.CustomCostFcn is OK. Analysis of user-provided model, cost, and constraint functions complete.

Analysis of Optimization Results

Find the optimal trajectories for the manipulated variables such that production ofCis maximized at the end of the batch process. To do so, use thenlmpcmovefunction.

fprintf('\nOptimization started...\n'); [~,~,Info] = nlmpcmove(nlmpcobj_Plan,x0,u0(1:2),zeros(1,3),u0(3)); fprintf(' Expected production of C (y1) is %g moles.\n',Info.Yopt(end,1)); fprintf(' First order optimality is satisfied (Info.ExitFlag = %i).\n',...Info.ExitFlag); fprintf('Optimization finished...\n');

Optimization started... Slack variable unused or zero-weighted in your custom cost function. All constraints will be hard. Expected production of C (y1) is 2.02353 moles. First order optimality is satisfied (Info.ExitFlag = 1). Optimization finished...

The discretized model uses a simple Euler integration, which could be inaccurate. To check this, integrate the model using theode15scommand for the calculated optimal MV trajectory.

Nstep = size(Info.Xopt,1) - 1; t = 0; X = x0'; t0 = 0;fori = 1:Nstep u_in = [Info.MVopt(i,1:2)'; c_Bin]; ODEFUN = @(t,x) fedbatch_StateFcn(x, u_in); TSPAN = [t0, t0+Ts]; Y0 = X(end,:)'; [TOUT,YOUT] = ode15s(ODEFUN,TSPAN,Y0); t = [t; TOUT(2:end)]; X = [X; YOUT(2:end,:)]; t0 = t0 + Ts;endnx = size(X,1); Y = zeros(nx,3);fori = 1:nx Y(i,:) = fedbatch_OutputFcn(X(i,:)',u_in)';endfprintf('\n Actual Production of C (y1) is %g moles.\n',X(end,2)-X(end,1)); fprintf(' Heat removal rate (y2) satisfies the upper bound.\n');

Actual Production of C (y1) is 2.0228 moles. Heat removal rate (y2) satisfies the upper bound.

In the top plot of the following figure, the actual production ofCagrees with the expected production ofCcalculated fromnlmpcmove. In the bottom plot, the heat removal rate never exceeds its hard constraint.

figure subplot(2,1,1) plot(t,Y(:,1),(0:Nstep)*Ts, Info.Yopt(:,1),'*') axis([0 0.5 0 Y(end,1) + 0.1]) legend({'Actual','Expected'},'location','northwest') title('Mol C in reactor (y1)') subplot(2,1,2) tTs = (0:Nstep)*Ts; t(end) = 0.5; plot(t,Y(:,2),“- - -”,[0 tTs(end)],1.5e5*ones(1,2),'r--') axis([0 0.5 0.8e5, 1.6e5]) legend({'q_r','Upper Bound'},'location','southwest') title('Heat removal rate (y2)')

Close examination of the heat removal rate shows that it can exhibit peaks and valleys between the sampling instants as reactant compositions change. Consequently, the heat removal rate exceeds the specified maximum of1.45e5(约0.35 t = h)但保持低于真实的格言um of1.5e5.

The following figure shows the optimal trajectory of planned adjustments in theBfeed rate (u1), and the reactor temperature (x4) and its setpoint (u2).

figure subplot(2,1,1) stairs(tTs,Info.MVopt(:,1)) title('Feed rate of B (u1)') subplot(2,1,2) plot(tTs,Info.MVopt(:,2),'*',t,X(:,4)-273.15,“- - -”,...[0 0.5],[20 20],'r--',[0 0.5],[50 50],'r--') axis([0 0.5 15 55]) title('Reactor temperature and its setpoint') legend({'Setpoint','Actual'},'location','southeast')

The trajectory begins with a relatively high feed rate, which increasesc_Band the resultingCproduction rate. To prevent exceeding the heat removal rate constraint, reactor temperature and feed rate must decrease. The temperature eventually hits its lower bound and stays there until the reactor is nearly full and theBfeed rate must go to zero. The temperature then increases to its maximum (to increaseCproduction) and finally drops slightly (to reduce D production, which is favored at higher temperatures).

The top plot of the following figure shows the consumption ofc_A, which tends to reduceCproduction. To compensate, the plan first increasesc_B, and when that is no longer possible (the reactor liquid volume must not exceed 1.1), the plan makes optimal use of the temperature. In the bottom plot of the following figure, the liquid volume never exceeds its upper bound.

figure subplot(2,1,1) c_A = X(:,1)./X(:,3); c_B = (c_Bin*X(:,3) + X(:,1) + V0*(c_B0 - c_A0 - c_Bin))./X(:,3); plot(t,[c_A, c_B]) legend({'c_A','c_B'},'location','west') subplot(2,1,2) plot(tTs,Info.Yopt(:,3)) title('Liquid volume')

Nonlinear MPC Design for Tracking the Optimal C Product Trajectory

To track the optimal trajectory of productCcalculated above, you design another nonlinear MPC controller with the same prediction model and constraints. However, use the standard quadratic cost and default horizons for tracking purposes.

To simplify the control task, assume that the optimal trajectory of theBfeed rate is implemented in the plant and the tracking controller considers it to be a measured disturbance. Therefore, the controller uses the reactor temperature setpoint as its only manipulated variable to track the desiredy1profile.

Create the tracking controller.

nlmpcobj_Tracking = nlmpc(4,3,'MV',2,'MD',[1,3]); nlmpcobj_Tracking.Ts = Ts; nlmpcobj_Tracking.Model = nlmpcobj_Plan.Model; nlmpcobj_Tracking.MV = nlmpcobj_Plan.MV(2); nlmpcobj_Tracking.OV = nlmpcobj_Plan.OV; nlmpcobj_Tracking.Weights.OutputVariables = [1 0 0];% track y1 onlynlmpcobj_Tracking.Weights.ManipulatedVariablesRate = 1e-6;% agressive MV

Zero weights are applied to one or more OVs because there are fewer MVs than OVs.

Obtain theCproduction (y1) reference signal from the optimal plan trajectory.

Cref = Info.Yopt(:,1);

Obtain the feed rate ofB(u1) from the optimal plan trajectory. The feed concentration ofB(u3) is a constant.

MD = [Info.MVopt(:,1) c_Bin*ones(N+1,1)];

First, run the tracking controller in nonlinear MPC mode.

[X1,Y1,MV1,et1] = fedbatch_Track(nlmpcobj_Tracking,x0,u0(2),N,Cref,MD); fprintf('\nNonlinear MPC: Elapsed time = %g sec. Production of C = %g mol\n',et1,Y1(end,1));

Nonlinear MPC: Elapsed time = 5.98252 sec. Production of C = 2.01754 mol

Second, run the controller as an adaptive MPC controller.

nlmpcobj_Tracking.Optimization.RunAsLinearMPC ='Adaptive'; [X2,Y2,MV2,et2] = fedbatch_Track(nlmpcobj_Tracking,x0,u0(2),N,Cref,MD); fprintf('\nAdaptive MPC: Elapsed time = %g sec. Production of C = %g mol\n',et2,Y2(end,1));

Adaptive MPC: Elapsed time = 3.0612 sec. Production of C = 2.01567 mol

Third, run the controller as a time-varying MPC controller.

nlmpcobj_Tracking.Optimization.RunAsLinearMPC ='TimeVarying'; [X3,Y3,MV3,et3] = fedbatch_Track(nlmpcobj_Tracking,x0,u0(2),N,Cref,MD); fprintf('\nTime-varying MPC: Elapsed time = %g sec. Production of C = %g mol\n',et3,Y3(end,1));

Time-varying MPC: Elapsed time = 2.50141 sec. Production of C = 2.02346 mol

在大多数MPC应用,线性MPClutions, such as Adaptive MPC and Time-varying MPC, provide performance that is comparable to the nonlinear MPC solution, while consuming less resources and executing faster. In these cases, nonlinear MPC often represents the best control results that MPC can achieve. By running a nonlinear MPC controller as a linear MPC controller, you can assess whether implementing a linear MPC solution is good enough in practice.

In this example, all three methods come close to the optimalCproduction obtained in the planning stage.

figure plot(Ts*(0:N),[Y1(:,1) Y2(:,1) Y3(:,1)]) title('Production of C') legend({'NLMPC','Adaptive','TimeVarying'},'location','northwest')

The unexpected result is that time-varying MPC produces moreCthan nonlinear MPC. The explanation is that the model linearization approaches used in the adaptive and time-varying modes result in a violation of the heat removal constraint, which results in a higherCproduction.

figure plot(Ts*(0:N),[Y1(:,2) Y2(:,2) Y3(:,2) 1.5e5*ones(N+1,1)]) title('Heat removal rate') legend({'NLMPC','Adaptive','TimeVarying','Constraint'},'location','southwest')

The adaptive MPC mode uses the plant states and inputs at the beginning of each control interval to obtain a single linear prediction model. This approach does not account for the known future changes in the feed rate, for example.

The time-varying method avoids this issue. However, at the start of the batch it assumes (by default) that the states will remain constant over the horizon. It corrects for this once it obtains its first solution (using data in theoptsvariable), but its initial choice of reactor temperature is too high, resulting in an earlyq_rconstraint violation.

References

[1] Srinivasan, B., S. Palanki, and D. Bonvin, "Dynamic optimization of batch processes I. Characterization of the nominal solution",Computers and Chemical Engineering, vol. 27 (2003), pp. 1-26.