function [] = forcesmpc_miso()
%% Control of a Multi-Input Single-Output Plant
% This example shows how to design model predictive controller with one
% measured output, one manipulated variable, one measured disturbance, and
% one unmeasured disturbance in a typical workflow.
% It also shows how to solve the MPC problems with unmeasured and measured
% disturbances using FORCESPRO.
%
% The example can be run with two different algorithms: The QP "general"
% algorithm and the "QP fast" algorithm.
%
% Guide for "QP general" algorithm:
% 1) run 'forcesMpc_miso()' in your Matlab terminal.
%
% Guide for "QP fast" algorithm:
% The "QP fast" algorithm tunes certain algorithm parameters according to
% the problem formulation for which it is designed. To do this one has to
% first collect data for tuning the solver. Collecting data is done using
% the following 3 commands:
% 1) run 'forcesMpcCache on' in your Matlab terminal.
% 2) run 'forcesMpc_miso()' in your Matlab terminal.
% 3) run 'forcesMpcCache off' in your Matlab terminal.
% After data has been collected you can run the "QP fast" algorithm as
% follows:
% 1) run 'forcesMpc_miso()' in your Matlab terminal.
%
% This example requires Model Predictive Control Toolbox from MathWorks.
% Copyright 2019-2023 The MathWorks, Inc. and Embotech AG, Zurich, Switzerland
clc;
close all;
%% Define Plant Model
% The discrete-time linear open-loop dynamic model is defined below with
% sample time |Ts|.
sys = ss(tf({1,1,1},{[1 .5 1],[1 1],[.7 .5 1]}));
Ts = 0.2;
model = c2d(sys,Ts);
%% Design MPC Controller
% Define type of input signals: the first signal is a manipulated variable,
% the second signal is a measured disturbance, the third one is an
% unmeasured disturbance.
model = setmpcsignals(model,'MV',1,'MD',2,'UD',3);
%%
% Create the controller object with sampling period, prediction and control
% horizons.
mpcobj = mpc(model,Ts,10,3);
%%
% Define constraints on the manipulated variable.
mpcobj.MV = struct('Min',-0.8,'Max',0.8,'RateMin',-10,'RateMax',10);
%%
% For unmeasured input disturbances, its model is an integrator driven by
% white noise with variance = 1000.
mpcobj.Model.Disturbance = tf(sqrt(1000),[1 0]);
%% Simulate Controller Using |mpcmove|
% First, obtain the discrete-time state-space matrices of the plant.
[A, B, C, D] = ssdata(model);
Tstop = 30; % Simulation time
%%
% Create options structure for FORCESPRO solver using the Sparse QP
% formulation.
options = mpcToForcesOptions('dense');
% Generate FORCESPRO solver
[coredata, statedata, onlinedata] = mpcToForces(mpcobj, options);
%%
% Store the closed-loop MPC trajectories in arrays |YY|, |UU|, and |XX|.
YY=[];
UU=[];
XX=[];
%% Simulate model with measured and unmeasured disturbance
simuSteps = round(Tstop/Ts);
% Reference signal
yref = 0.5*ones(simuSteps, 1);
% Measured disturbance signal
mdis = zeros(simuSteps, 1);
mdis(50:end) = 1;
% Initial state of the plant
x = [0 0 0 0 0]';
for t=1:simuSteps
% Get online preview signals
onlinedata.signals.ref = yref(t:min(t+mpcobj.p-1,simuSteps));
onlinedata.signals.md = mdis(t:min(t+mpcobj.p,simuSteps));
% Store states
XX = [XX,x];
% Define unmeasured disturbance signal
d = 0;
if t*Ts>=20
d = -0.5;
end
% Plant equations: output update (no feedthrough from MV to Y)
y = C*x + D(:,2)*mdis(t) + D(:,3)*d;
YY = [YY,y];
% Compute MPC action using FORCESPRO
onlinedata.signals.ym = y;
[mv, statedata, ~] = mpcmoveForces(coredata, statedata, onlinedata);
% Plant equations: state update
x = A*x + B*[mv;mdis(t);d];
UU = [UU,mv];
end
%% Plot results
figure
subplot(2,1,1);
plot(Ts:Ts:Tstop,YY,Ts:Ts:Tstop,yref);
grid on; xlabel('Time (sec)'); title('Output');
subplot(2,1,2);
plot(0:Ts:Tstop-Ts,UU,0:Ts:Tstop-Ts,-0.8*ones(1,simuSteps),0:Ts:Tstop-Ts,0.8*ones(1,simuSteps));
grid on; xlabel('Time (sec)'); title('Input');
end