Model Predictive Control of Multi

您所在的位置：网站首页 › plantmodel › Model Predictive Control of Multi

Model Predictive Control of Multi

#Model Predictive Control of Multi| 来源: 网络整理| 查看: 265

Define Plant Model

Define a plant model. For this example, use continuous-time transfer functions from each input to the output.

plantTF = tf({1,1,1},{[1 .5 1],[1 1],[.7 .5 1]}) % define and display tf object plantTF = From input 1 to output: 1 --------------- s^2 + 0.5 s + 1 From input 2 to output: 1 ----- s + 1 From input 3 to output: 1 ------------------- 0.7 s^2 + 0.5 s + 1 Continuous-time transfer function.

For this example, explicitly convert the plant to a discrete-time state-space form before passing it to the MPC controller creation function.

The controller creation function can accept either continuous-time or discrete-time plants. During initialization, a continuous-time plant (in any format) is automatically converted into a discrete-time state-space model using the zero-order hold (ZOH) method. Delays, if present, are incorporated in the state-space model.

You can convert the plant to discrete-time yourself when you need the discrete-time system matrices for analysis or simulation (as in this example) or when you want to use a discrete-time conversion method other than ZOH.

plantCSS = ss(plantTF); % transfer function to continuous state space Ts = 0.2; % specify a sample time of 0.2 seconds plantDSS = c2d(plantCSS,Ts) % convert to discrete-time state space, using ZOH plantDSS = A = x1 x2 x3 x4 x5 x1 0.8862 -0.1891 0 0 0 x2 0.1891 0.9807 0 0 0 x3 0 0 0.8187 0 0 x4 0 0 0 0.841 -0.2637 x5 0 0 0 0.1846 0.9729 B = u1 u2 u3 x1 0.1891 0 0 x2 0.01929 0 0 x3 0 0.1813 0 x4 0 0 0.1846 x5 0 0 0.01899 C = x1 x2 x3 x4 x5 y1 0 1 1 0 1.429 D = u1 u2 u3 y1 0 0 0 Sample time: 0.2 seconds Discrete-time state-space model.

By default, the software assumes that all the plant input signals are manipulated variables. To specify the signal types, such as measured and unmeasured disturbances, use the setmpcsignals function. In this example, the first input signal is a manipulated variable, the second is a measured disturbance, and the third is an unmeasured disturbance. This information is stored in the plant model plantDSS and later used by the MPC controller.

plantDSS = setmpcsignals(plantDSS,'MV',1,'MD',2,'UD',3); % specify signal types Design MPC Controller

Create the controller object, specifying the sample time, as well as the prediction and control horizons (10 and 3 steps, respectively).

mpcobj = mpc(plantDSS,Ts,10,3); -->"Weights.ManipulatedVariables" is empty. Assuming default 0.00000. -->"Weights.ManipulatedVariablesRate" is empty. Assuming default 0.10000. -->"Weights.OutputVariables" is empty. Assuming default 1.00000.

Since you have not specified the weights of the quadratic cost function to be minimized, the controller uses their default values (0 for manipulated variables, 0.1 for manipulated variable rates, and 1 for the output variables). Also, at this point the MPC problem is still unconstrained as you have not specified any constraint yet.

Define hard constraints on the manipulated variable.

mpcobj.MV = struct('Min',0,'Max',1,'RateMin',-10,'RateMax',10);

You can use input and output disturbance models to define the dynamic characteristics of additive input and output unmeasured disturbances. These models allow the controller to better reject such disturbances, if they occur at run time. By default, to be able to better reject step-like disturbances, mpc uses an integrator as disturbance model for:

Each unmeasured disturbance input and

Each unmeasured disturbance acting on each measured outputs

unless doing so causes a violation of state observability.

The MPC object also has a noise model that specifies the characteristics of the additive noise that is expected on the measured output variables. By default, this, is assumed to be a unit static gain, which is equivalent to assuming that the controller expects the measured output variables to be affected, at run time, by white noise (with a covariance matrix that depends on the input matrices of the whole prediction model). For more information, see MPC Prediction Models.

Display the input disturbance model. As expected, it is a discrete-time integrator.

getindist(mpcobj) -->The "Model.Disturbance" property is empty: Assuming unmeasured input disturbance #3 is integrated white noise. Assuming no disturbance added to measured output #1. -->"Model.Noise" is empty. Assuming white noise on each measured output. ans = A = x1 x1 1 B = UD1-wn x1 0.2 C = x1 UD1 1 D = UD1-wn UD1 0 Sample time: 0.2 seconds Discrete-time state-space model.

Display the output disturbance model.

getoutdist(mpcobj) ans = Empty state-space model.

Specify the disturbance model for the unmeasured input as an integrator driven by white noise with variance 1000.

mpcobj.Model.Disturbance = tf(sqrt(1000),[1 0]);

Display the input disturbance model again to verify that it changed.

getindist(mpcobj) ans = A = x1 x1 1 B = Noise#1 x1 0.8 C = x1 UD1 7.906 D = Noise#1 UD1 0 Sample time: 0.2 seconds Discrete-time state-space model.

Display the MPC controller object mpcobj to review its properties.

mpcobj MPC object (created on 19-Aug-2023 23:11:21): --------------------------------------------- Sampling time: 0.2 (seconds) Prediction Horizon: 10 Control Horizon: 3 Plant Model: -------------- 1 manipulated variable(s) -->| 5 states | | |--> 1 measured output(s) 1 measured disturbance(s) -->| 3 inputs | | |--> 0 unmeasured output(s) 1 unmeasured disturbance(s) -->| 1 outputs | -------------- Indices: (input vector) Manipulated variables: [1 ] Measured disturbances: [2 ] Unmeasured disturbances: [3 ] (output vector) Measured outputs: [1 ] Disturbance and Noise Models: Output disturbance model: default (type "getoutdist(mpcobj)" for details) Input disturbance model: user specified (type "getindist(mpcobj)" for more details) Measurement noise model: default (unity gain after scaling) Weights: ManipulatedVariables: 0 ManipulatedVariablesRate: 0.1000 OutputVariables: 1 ECR: 100000 State Estimation: Default Kalman Filter (type "getEstimator(mpcobj)" for details) Constraints: 0 "Model.Noise" is empty. Assuming white noise on each measured output. -->In closed-loop simulation, the "MVSignal" property of "mpcsimopt" object is ignored. -->Converting model to discrete time.

Plot the plant output response as well as the plant states estimated by the new observer.

figure; % create figure subplot(2,1,1) % create upper subplot axis plot(t,y) % plot y versus time title('Plant Output'); % add title to upper plot ylabel('y') % add a label to the upper y axis grid % add grid to upper plot subplot(2,1,2) % create lower subplot axis plot(t,xc.Plant) % plot xc.Plant versus time title('Estimated Plant States'); % add title to lower plot xlabel('Time (seconds)') % add a label to the lower x axis ylabel('xc') % add a label to the lower y axis grid % add grid to lower plot

As expected, the controller states are different from the ones previously plotted, and the overall closed-loop response is somewhat slower.

Simulate Controller in Closed Loop Using mpcmove

For more general applications, you can simulate an MPC controller in a for loop using the mpcmove function. Using this function, you can run simulations with the following features.

Nonlinear or time-varying plants

Constraints or weights that vary at run time

Disturbance or reference signals that are not known before running the simulation

If your plant is continuous, you can either convert it to discrete time before simulating or you can use a numerical integration algorithm (such as forward Euler or ode45) to simulate it in a closed loop using mpcmove. For example, you can calculate the plant state at the next control interval using the following methods:

Discrete time plant x(t+1)=f(x(t),u(t)): x = f(x,u), (typically x = A*x + B*u for linear plant models)

Continuous time plant dx(t)/dt=f(x(t),u(t)), sample time Ts, Euler method: x = x + f(x,u)*Ts

Continuous time plant as above, using ode45: [~,xhist] = ode45(@(t,xode) f(xode,u),[0 Ts],x); x = xhist(end);

In the third case, ode45 starts from the initial condition x and simulates the plant for Ts seconds, under a constant control signal u. The last value of the resulting state signal xhist is the plant state at the next control interval.

First, obtain the discrete-time state-space matrices of the plant, and define the simulation time and initial states for plant and controller.

[A,B,C,D] = ssdata(plantDSS); % discrete-time plant plant ss matrices Tstop = 30; % simulation time x = [0 0 0 0 0]'; % initial state of the plant xmpc = mpcstate(mpcobj); % get handle to controller state r = 1; % output reference signal

Display the initial state of the controller. The state, which is an mpcstate object, contains the controller states only at the current time. Specifically: * |xmpc.Plant| is the current value of the estimated plant states. * |xmpc.Disturbance| is the current value of the disturbance models states. * |xmpc.Noise| is the current value of the noise models states. * |xmpc.LastMove| is the last value of the manipulated variable. * |xmpc.Covariance| is the current value of the estimator covariance matrix.

xmpc % display controller states MPCSTATE object with fields Plant: [0 0 0 0 0] Disturbance: 0 Noise: [1x0 double] LastMove: 0 Covariance: []

Note that xmpc is a handle object, which always points to the current state of the controller. Since mpcmove updates the internal plant state when a new control move is calculated, you do not need to update xmpc, which always points to the current (hence updated) state.

isa(xmpc,'handle') ans = logical 1

Define workspace arrays YY and UU to store output and input signals, respectively, so that you can plot them after the simulation.

YY=[]; UU=[];

Run the simulation loop.

for k=0:round(Tstop/Ts)-1 % Define measured disturbance signal v(k). You can specify a more % complex dependence on time or previous states here, if needed. v = 0; if k*Ts>=10 % raising to 1 after 10 seconds v = 1; end % Define unmeasured disturbance signal d(k). You can specify a more % complex dependence on time or previous states here, if needed. d = 0; if k*Ts>=20 % falling to -0.5 after 20 seconds d = -0.5; end % Plant equations: current output % If you have a more complex plant output behavior (including, for example, % model mismatch or nonlinearities) you can to simulate it here. % Note that there cannot be any direct feedthrough between u and y. y = C*x + D(:,2)*v + D(:,3)*d; % calculate current output (D(:,1)=0) YY = [YY,y]; % store current output % Note, if the plant had a non-zero operating point the output would be: % y = mpcobj.Model.Nominal.Y + C*(x-mpcobj.Model.Nominal.X) + D(:,2)*v + D(:,3)*d; % Compute the MPC action (u) and update the internal controller states. % Note that you do not need the update xmpc because it always points % to the current controller states. u = mpcmove(mpcobj,xmpc,y,r,v); % xmpc,y,r,v are values at current step k UU = [UU,u]; % store current input % Plant equations: state update % You can simulate a more complex plant state behavior here, if needed. x = A*x + B(:,1)*u + B(:,2)*v + B(:,3)*d; % update state % Note, if the plant had a non-zero operating point the state update would be: % x = mpcobj.Model.Nominal.X + mpcobj.Model.Nominal.DX + ... % A*(x-mpcobj.Model.Nominal.X) + B(:,1)*(u-mpcobj.Model.Nominal.U(1)) + ... % B(:,2)*v + B(:,3)*d; end

Plot the results.

figure % create figure subplot(2,1,1) % create upper subplot axis plot(0:Ts:Tstop-Ts,YY) % plot YY versus time ylabel('y') % add a label to the upper y axis grid % add grid to upper plot title('Output') % add title to upper plot subplot(2,1,2) % create lower subplot axis plot(0:Ts:Tstop-Ts,UU) % plot UU versus time ylabel('u') % add a label to the lower y axis xlabel('Time (seconds)') % add a label to the lower x axis grid % add grid to lower plot title('Input') % add title to lower plot

To check the optimal predicted trajectories at any point during the simulation, you can use the second output argument of mpcmove. For this example, assume you start from the current state (x and xmpc). Also assume that, from this point until the end of the horizon, the reference set-point is 0.5 and the disturbance is 0. Simulate the controller and return the info structure.

r = 0.5; % reference v = 0; % disturbance [~,info] = mpcmove(mpcobj,xmpc,y,r,v); % solve over prediction horizon

Display the info variable.

info info = struct with fields: Uopt: [11x1 double] Yopt: [11x1 double] Xopt: [11x6 double] Topt: [11x1 double] Slack: 0 Iterations: 1 QPCode: 'feasible' Cost: 0.1399

info is a structure containing the predicted optimal sequences of manipulated variables, plant states, and outputs over the prediction horizon. mpcmove calculated this sequence, together with the optimal first move, by solving a quadratic optimization problem to minimize the cost function. The plant states and outputs in info result from applying the optimal manipulated variable sequence directly to mpcobj.Model.Plant, in an open-loop fashion. Due to the presence of noise, unmeasured disturbances, and uncertainties, this open-loop optimization process is not equivalent to simulating the closed loop consisting of the plant, estimator and controller using either the sim command or mpcmove iteratively in a for loop.

Extract the predicted optimal trajectories.

topt = info.Topt; % time yopt = info.Yopt; % predicted optimal plant model outputs uopt = info.Uopt; % predicted optimal mv sequence

Since the optimal sequence values are constant across each control step, plot the trajectories using a stairstep plot.

figure % create new figure subplot(2,1,1) % create upper subplot stairs(topt,yopt) % plot yopt in a stairstep graph title('Optimal Sequence of Predicted Outputs') % add title to upper subplot grid % add grid to upper subplot subplot(2,1,2) % create lower subplot stairs(topt,uopt) % plot uopt in a stairstep graph title('Optimal Sequence of Manipulated Variables') % add title to upper subplot grid % add grid to upper subplot

Linear Representation of MPC Controller

When the constraints are not active, the MPC controller behaves like a linear controller. Note that for a finite-time unconstrained Linear Quadratic Regulator problem with a finite non-receding horizon, the value function is time-dependent, which causes the optimal feedback gain to be time varying. In contrast, in MPC the horizon has a constant length because it is always receding, resulting in a time-invariant value function and consequently a time-invariant optimal feedback gain.

You can get the state-space form of the MPC controller.

LTI = ss(mpcobj,'rv'); % get state-space representation

Get the state-space matrices to simulate the linearized controller.

[AL,BL,CL,DL] = ssdata(LTI); % get state-space matrices

Initialize variables for a closed-loop simulation of both the original MPC controller without constraints and the linearized controller.

mpcobj.MV = []; % remove input constraints mpcobj.OV = []; % remove output constraints Tstop = 5; % set simulation time y = 0; % set nitial measured output r = 1; % set output reference set-point (constant) u = 0; % set previous (initial) input command x = [0 0 0 0 0]'; % set initial state of plant xmpc = mpcstate(mpcobj); % set initial state of unconstrained MPC controller xL = zeros(size(BL,1),1); % set initial state of linearized MPC controller YY = []; % define workspace array to store plant outputs Assuming no disturbance added to measured output #1. -->"Model.Noise" is empty. Assuming white noise on each measured output.

Simulate both controllers in a closed loop with the same plant model.

for k = 0:round(Tstop/Ts)-1 YY = [YY,y]; % store current output for plotting purposes % Define measured disturbance signal v(k). v = 0; if k*Ts>=10 v = 1; % raising to 1 after 10 seconds end % Define unmeasured disturbance signal d(k). d = 0; if k*Ts>=20 d = -0.5; % falling to -0.5 after 20 seconds end % Compute the control actions of both (unconstrained) MPC and linearized MPC uMPC = mpcmove(mpcobj,xmpc,y,r,v); % unconstrained MPC (also updates xmpc) u = CL*xL + DL*[y;r;v]; % unconstrained linearized MPC % Compare the two control actions dispStr(k+1) = {sprintf(['t=%5.2f, u=%7.4f (provided by LTI), u=%7.4f' ... ' (provided by MPCOBJ)'],k*Ts,u,uMPC)}; %#ok % Update state of unconstrained linearized MPC controller xL = AL*xL + BL*[y;r;v]; % Update plant state x = A*x + B(:,1)*u + B(:,2)*v + B(:,3)*d; % Calculate plant output y = C*x + D(:,1)*u + D(:,2)*v + D(:,3)*d; % D(:,1)=0 end

Display the character arrays containing the control actions.

for k=0:round(Tstop/Ts)-1 disp(dispStr{k+1}); % display each string as k increases end t= 0.00, u= 5.2478 (provided by LTI), u= 5.2478 (provided by MPCOBJ) t= 0.20, u= 3.0134 (provided by LTI), u= 3.0134 (provided by MPCOBJ) t= 0.40, u= 0.2281 (provided by LTI), u= 0.2281 (provided by MPCOBJ) t= 0.60, u=-0.9952 (provided by LTI), u=-0.9952 (provided by MPCOBJ) t= 0.80, u=-0.8749 (provided by LTI), u=-0.8749 (provided by MPCOBJ) t= 1.00, u=-0.2022 (provided by LTI), u=-0.2022 (provided by MPCOBJ) t= 1.20, u= 0.4459 (provided by LTI), u= 0.4459 (provided by MPCOBJ) t= 1.40, u= 0.8489 (provided by LTI), u= 0.8489 (provided by MPCOBJ) t= 1.60, u= 1.0192 (provided by LTI), u= 1.0192 (provided by MPCOBJ) t= 1.80, u= 1.0511 (provided by LTI), u= 1.0511 (provided by MPCOBJ) t= 2.00, u= 1.0304 (provided by LTI), u= 1.0304 (provided by MPCOBJ) t= 2.20, u= 1.0053 (provided by LTI), u= 1.0053 (provided by MPCOBJ) t= 2.40, u= 0.9920 (provided by LTI), u= 0.9920 (provided by MPCOBJ) t= 2.60, u= 0.9896 (provided by LTI), u= 0.9896 (provided by MPCOBJ) t= 2.80, u= 0.9925 (provided by LTI), u= 0.9925 (provided by MPCOBJ) t= 3.00, u= 0.9964 (provided by LTI), u= 0.9964 (provided by MPCOBJ) t= 3.20, u= 0.9990 (provided by LTI), u= 0.9990 (provided by MPCOBJ) t= 3.40, u= 1.0002 (provided by LTI), u= 1.0002 (provided by MPCOBJ) t= 3.60, u= 1.0004 (provided by LTI), u= 1.0004 (provided by MPCOBJ) t= 3.80, u= 1.0003 (provided by LTI), u= 1.0003 (provided by MPCOBJ) t= 4.00, u= 1.0001 (provided by LTI), u= 1.0001 (provided by MPCOBJ) t= 4.20, u= 1.0000 (provided by LTI), u= 1.0000 (provided by MPCOBJ) t= 4.40, u= 0.9999 (provided by LTI), u= 0.9999 (provided by MPCOBJ) t= 4.60, u= 1.0000 (provided by LTI), u= 1.0000 (provided by MPCOBJ) t= 4.80, u= 1.0000 (provided by LTI), u= 1.0000 (provided by MPCOBJ)

Plot the results.

figure % create figure plot(0:Ts:Tstop-Ts,YY) % plot plant outputs grid % add grid title('Unconstrained MPC control: Plant output') % add title xlabel('Time (seconds)') % add label to x axis ylabel('y') % add label to y axis

Running a closed-loop simulation in which all controller constraints are turned off is easier using sim, as you just need to specify 'off' in the Constraint field of the related mpcsimopt simulation option object.

SimOptions = mpcsimopt; % create simulation options object SimOptions.Constraints = 'off'; % remove all MPC constraints SimOptions.UnmeasuredDisturbance = d; % unmeasured input disturbance sim(mpcobj,Nf,r,v,SimOptions); % run closed-loop simulation

Simulate Using Simulink

You can also simulate your MPC controller in Simulink.

To compare results, recreate the MPC object with the constraints you use in the Design MPC Controller section, and the default estimator.

mpcobj = mpc(plantDSS,Ts,10,3); mpcobj.MV = struct('Min',0,'Max',1,'RateMin',-10,'RateMax',10); mpcobj.Model.Disturbance = tf(sqrt(1000),[1 0]); -->"Weights.ManipulatedVariables" is empty. Assuming default 0.00000. -->"Weights.ManipulatedVariablesRate" is empty. Assuming default 0.10000. -->"Weights.OutputVariables" is empty. Assuming default 1.00000.

Obtain the state-space matrices of the continuous-time plant.

[A,B,C,D] = ssdata(plantCSS); % get state-space realization

Open the mpc_miso Simulink model for closed-loop simulation. The plant model is implemented with a continuous state-space block.

open_system('mpc_miso')

The plant input signals u(t), v(t), and d(t) represent the manipulated variable, measured input disturbance, and unmeasured input disturbance, respectively, while y(t) is the measured output. The block parameters are the matrices forming the state-space realization of the continuous-time plant, and the initial conditions for the five states. The MPC controller is implemented with an MPC Controller block, which has the workspace MPC object mpcobj as a parameter, the manipulated variable as the output, and the measured plant output, reference signal, and measured plant input disturbance, respectively, as inputs. The four Scope blocks plot the five loop signals, which are also saved (except for the reference signal) by four To-Workspace blocks.

Simulate the closed loop system using the simulink sim command. Note that this command (which simulates a Simulink model, and is equivalent to clicking the "Run" button in the model) is different from the sim command provided by the MPC toolbox (which instead simulates an MPC controller in a loop with an LTI plant).

sim('mpc_miso') Assuming no disturbance added to measured output #1. -->"Model.Noise" is empty. Assuming white noise on each measured output.

To show the simulation results, open the four Scope windows

open_system('mpc_miso/MV') open_system('mpc_miso/Outputs//References') open_system('mpc_miso/MD') open_system('mpc_miso/UD')

The plots in the Scope windows are equivalent to the ones in the Simulate Closed-Loop Response Using the sim Command and Simulate Closed-Loop Response with Model Mismatch sections, with minor differences due to the fact that in Simulate Closed-Loop Response Using the sim Command the unmeasured disturbance signal is zero, and that in Simulate Closed-Loop Response with Model Mismatch you add noise to the plant input and output. Also note that, while the MPC sim command internally discretizes any continuous plant model using the ZOH method, Simulink typically uses an integration algorithm (in this example ode45) to simulate the closed loop when a continuous-time block is present.

Run Simulation with Sinusoidal Output Noise.

Assume output measurements are affected by a sinusoidal disturbance (a single tone sensor noise) of frequency 0.1 Hz.

omega = 2*pi/10; % disturbance radial frequency

Open the mpc_misonoise Simulink model, which is similar to the mpc_miso model except for the sinusoidal disturbance added to the measured output. Also, the simulation time is longer and the unmeasured disturbance begins before the measured disturbance.

open_system('mpc_misonoise') % open new Simulink model

Since this noise is expected, you can specify a noise model to help the state estimator ignore it. Doing so improves the disturbance rejection capabilities of the controller.

mpcobj.Model.Noise = 0.5*tf(omega^2,[1 0 omega^2]); % measurement noise model

Revise the MPC design by specifying a disturbance model on the unmeasured input as a white Gaussian noise with zero mean and variance 0.1.

setindist(mpcobj,tf(0.1)); % static gain

In this case, you cannot have integrators as disturbance model on both the unmeasured input and the output, because this violates state observability. Therefore when you specify a static gain for the input disturbance model, an output disturbance model consisting in a discretized integrator is automatically added to the controller. This output disturbance model helps the controller to reject step-like and slowly varying disturbances at the output.

getoutdist(mpcobj) -->Assuming output disturbance added to measured output #1 is integrated white noise. -->A feedthrough channel in NoiseModel was inserted to prevent problems with estimator design. ans = A = x1 x1 1 B = u1 x1 0.2 C = x1 MO1 1 D = u1 MO1 0 Sample time: 0.2 seconds Discrete-time state-space model.

Large measurement noise can decrease the accuracy of the state estimates. To make the controller less aggressive, and decrease its noise sensitivity, decrease the weight on the output variable tracking.

mpcobj.weights = struct('MV',0,'MVRate',0.1,'OV',0.005); % new weights

To give the Kalman filter more time to successfully estimate the states, increase the prediction horizon to 40.

mpcobj.predictionhorizon = 40; % new prediction horizon

Run the simulation for 145 seconds.

sim('mpc_misonoise',145) % the second argument is the simulation duration -->Assuming output disturbance added to measured output #1 is integrated white noise. -->A feedthrough channel in NoiseModel was inserted to prevent problems with estimator design.

To show the simulation results, open the four Scope windows

open_system('mpc_misonoise/MV') open_system('mpc_misonoise/Outputs//References//Noise') open_system('mpc_misonoise/MD') open_system('mpc_misonoise/UD')

The Kalman filter successfully learns to ignore the measurement noise after 50 seconds. The unmeasured and measured disturbances get rejected in a 10 to 20 second timespan. As expected, the manipulated variable stays in the interval between 0 and 1.

bdclose all % close all open Simulink models without saving any change close all % close all open figures

【本文地址】

公司简介

联系我们