Model Predictive Control

Suba Thomas & Sergio Vargas, Wolfram Research

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

Abstract

This talk will showcase the recently added functionality to design model predictive controllers. The formulation of the problem as a multiparametric optimization problem, along with workflows and related functionality will be presented.

Introduction

The control problem

Model

x(k+1)==a.x(k)+b.u(k)
y(k)==c.u(k)+d.x(k)

• StateSpaceModel
• AffineStateSpaceModel
• NonlinearStateSpaceModel
• SystemModel

Cost

(Subscript[J, 2])=(∑_(k=0)^(η-1))((Transpose[x(k)]).q.x(k)+(Transpose[u(k)]).r.u(k)+2(Transpose[x(k)]).p.u(k))+(Transpose[x(η)]).q.x(η)

• regulation or tracking
• 1-, squared 2-, and ∞ norms
• finite or infinite horizon

Constraints

xmin≤(Subscript[x, i])≤xmax
umin≤(Subscript[u, i])≤umax

In a nutshell

• model + cost + constraints
	• a optimization problem is solved at each instant
	• reformulate it as a parametric optimization problem to get a closed-			form solution

MPC as a parametric optimization problem

Parametric optimization

res=Minimize[(v1^2)+(v2^2)+(p^2),-1<=v1+v2<=1&&-5<=v1-p<=5&&-10<=p<=10,{v1,v2},Reals]//Expand

colorFn:=ColorData[97,Which[#1<=-6,1,#1<=-5,2,#1<5,3,#1<=6,4,True,5]]&

With[{plotOpts={Frame->True,Axes->False,PlotStyle->Thick,ColorFunction->colorFn,ColorFunctionScaling->False}},{Plot[v1/.res[[2]],{p,-10,10},plotOpts],Plot[v2/.res[[2]],{p,-10,10},plotOpts]}]

Plot[res[[1]],{p,-10,10},Frame->True,Axes->False,GridLines->Automatic,PlotStyle->Thick,ColorFunction->colorFn,ColorFunctionScaling->False]

MPC

The control problem.

model= x[k+1]==x[k]+u[k];

cost =(∑_(k=0)^2) (((x[k])^2)+((u[k])^2))

cons= Table[Splice[{-5<=x[k]<=5,-1<=u[k]<=1}],{k,0,2}]

It’s reformulation.

modelRule=x[k_?Positive]:> x[k-1]+u[k-1]

costNew=cost//.modelRule

consNew=cons//.modelRule

Variables[{costNew,consNew}]

res=Minimize[costNew,consNew, {u[0],u[1],u[2]}]//Simplify

modelSim=Table[x[i]//.modelRule,{i,0,2}]

Table[modelSim/.res[[2]]/.x[0]->RandomReal[{-5,5}],20];
ListLinePlot[%,PlotMarkers->Automatic,DataRange->{0,2},AxesLabel->{"k","x[k]"}]

MPC problems

Regulator problems

(Subscript[J, 2])=(∑_(k=0)^(η-1))((Transpose[(Transpose[x(k)]).q.x(k)+(Subscript[u, f])(k)]).r.(Subscript[u, f])(k)+2 (Transpose[x(k)]).p.(Subscript[u, f])(k))+(Transpose[x(η)]).q.x(η)

(Subscript[J, 1])=(∑_(k=0)^(η-1))(||q.x(k)(Subscript[||, 1])+||r.(Subscript[u, f])(k)(Subscript[||, 1])+2||(Subscript[p, x]).x(k)+(Subscript[p, u]).(Subscript[u, f])(k)(Subscript[||, 1]))+||q.x(η)(Subscript[||, 1])

(Subscript[J, ∞])=(∑_(k=0)^(η-1))(||q.x(k)(Subscript[||, ∞])+||r.(Subscript[u, f])(k)(Subscript[||, ∞])+2||(Subscript[p, x]).x(k)+(Subscript[p, u]).(Subscript[u, f])(k)(Subscript[||, ∞]))+||q.x(η)(Subscript[||, ∞])

Tracking problems

(Subscript[J, 2])=(∑_(k=0)^(η-1))((Transpose[e(k)]).q.e(k)+(Transpose[(Subscript[Δu, f])(k)]).r.(Subscript[Δu, f])(k)+2 (Transpose[e(k)]).p.(Subscript[Δu, f])(k))+(Transpose[e(η)]).q.e(η)

(Subscript[J, 1])=(∑_(k=0)^(η-1))(||q.e(k)(Subscript[||, 1])+||r.(Subscript[Δu, f])(k)(Subscript[||, 1])+2||(Subscript[p, e]).e(k)+(Subscript[p, u]).(Subscript[Δu, f])(k)(Subscript[||, 1]))+||q.e(η)(Subscript[||, 1])

(Subscript[J, ∞])=(∑_(k=0)^(η-1))(||q.e(k)(Subscript[||, ∞])+||r.(Subscript[Δu, f])(k)(Subscript[||, ∞])+2||(Subscript[p, e]).e(k)+(Subscript[p, u]).(Subscript[Δu, f])(k)(Subscript[||, ∞]))+||q.e(η)(Subscript[||, ∞])

Infinite horizon problems

(Subscript[J, 2])=(∑_(k=0)^∞)((Transpose[(Transpose[x(k)]).q.x(k)+(Subscript[u, f])(k)]).r.(Subscript[u, f])(k)+2 (Transpose[x(k)]).p.(Subscript[u, f])(k))

ModelPredictiveController

ModelPredictiveController[sspec,cost,cons]

System specifications

({{"InputModel"}, {"FeedbackInputs"}, {"TrackedOutputs"}})

Cost

({{"Regulator", "Tracker"}, {"StateWeight", "ErrorWeight"}, {"InputWeight", "InputIncrementWeight"}, {"StateInputWeight", "ErrorInputIncrementWeight"}, {"Horizon", "Horizon"}, {"Norm", "Norm"}})

Constraints

xmin≤(Subscript[x, i])≤xmax
min≤β+(Subscript[α, n]).(Subscript[x, n])+(Subscript[α, 1]).(Subscript[x, 1])+…≤max
xmin≤(Subscript[x, i])(k)≤xmax
Δxmin≤(Subscript[x, i])((Subscript[k, 2]))-(Subscript[x, i])((Subscript[k, 1]))≤Δxmax

Finite horizon 2 norm regulator

(({{(Subscript[x, 1])[k+1]}, {(Subscript[x, 2])[k+1]}}))==(({{1, 1}, {0, 1}})).(({{(Subscript[x, 1])[k]}, {(Subscript[x, 2])[k]}}))+(({{1}, {0.5}})).u(k)

ssm=StateSpaceModel[{(({{1, 1}, {0, 1}})),(({{1}, {0.5}}))},{(Subscript[x, 1]),(Subscript[x, 2])},u,SamplingPeriod->0.1]

(∑_(k=0)^3)({(Subscript[x, 1])(k),(Subscript[x, 2])(k)}.q.{(Subscript[x, 1])(k),(Subscript[x, 2])(k)}+{u(k)}.r.{u(k)})+{(Subscript[x, 1])(4),(Subscript[x, 2])(4)}.q.{(Subscript[x, 1])(4),(Subscript[x, 2])(4)}

{𝒩,q,r}= {4,DiagonalMatrix[{0.5,0}],{{0.5}}};

cost=<|"Horizon"->𝒩,"StateWeight"-> q,"InputWeight"-> r|>

cons=-5≤(Subscript[x, 1])≤5∧-5≤(Subscript[x, 2])≤5∧-1≤u≤1;

cd=ModelPredictiveController[ssm,cost,cons]

csys=cd["ClosedLoopSystem"]

OutputResponse[csys,PadRight[{0.1,0,-0.1,0,0,0.1},20]];
ListStepPlot[%,PlotRange->All,DataRange->{0, 0.1 20}]

Finite horizon 2 norm tracker

(({{(Subscript[x, 1])[k+1]}, {(Subscript[x, 2])[k+1]}}))==(({{1, 1}, {0, 1}})).(({{(Subscript[x, 1])[k]}, {(Subscript[x, 2])[k]}}))+(({{1}, {0.5}})).u(k)

y[k]==(({{1, 0}})).(({{(Subscript[x, 1])[k]}, {(Subscript[x, 2])[k]}}))

ssm=StateSpaceModel[{(({{1, 1}, {0, 1}})),(({{1}, {0.5}})),(({{1, 0}}))},{(Subscript[x, 1]),(Subscript[x, 2])},u,SamplingPeriod->0.1]

(∑_(k=0)^3)({(Subscript[e, 1])(k),(Subscript[e, 2])(k)}.qe.{(Subscript[e, 1])(k),(Subscript[e, 2])(k)}+{Δu(k)}.rΔ.{Δu(k)})+{(Subscript[e, 1])(4),(Subscript[e, 2])(4)}.qe.{(Subscript[e, 1])(4),(Subscript[e, 2])(4)}

{𝒩,qe,rΔ}= {4,{{10}},{{1}}};

cost=<|"Horizon"->𝒩,"ErrorWeight"-> qe,"InputIncrementWeight"->rΔ|>;

cons=-5≤(Subscript[x, 1])≤5∧-5≤(Subscript[x, 2])≤5∧-1≤u≤1;

cd=ModelPredictiveController[ssm,cost,cons]

csys=cd["ClosedLoopSystem"]

OutputResponse[csys,PadRight[{-1,2,3},30,2]];
ListStepPlot[%,PlotRange->All,DataRange->{0,0.1 30}]

ref=Table[Sin[t],{t,0,10,0.1}];
or=OutputResponse[csys,ref];

ListStepPlot[{ref,or[[1]]},PlotRange->All,DataRange->{0,0.1 20},PlotLegends->{"ref.", "act."}]

Infinite horizon 2 norm regulator

cd=ModelPredictiveController[(StateSpaceModelFull[1, 1, 0, 0, 1, 1, 0.1, ,{(Subscript[x, 1]),(Subscript[x, 2])},{u}, ]),<|"Horizon"->∞,"StateWeight"->IdentityMatrix[2],"InputWeight"->{{0.1}}|>,-3≤(Subscript[x, 1])≤3∧-3≤(Subscript[x, 2])≤3∧-1≤u≤1]

cd["FeedbackGainsModel"]//Simplify

reg=cd["FeedbackRegions"];
regPlot=RegionPlot[Evaluate[List@@reg[[1]]],{(Subscript[x, 1]),-3,3},{(Subscript[x, 2]),-3,3},PlotPoints->60,PlotTheme->"Web"]

csys=cd["ClosedLoopSystem"]

Manipulate[Show[regPlot,Graphics[{Thick,Arrow[OutputResponse[{csys,<|1->ic|>},
PadRight[{0},10]]\[Transpose]]}]],{ic,Locator},ImageMargins->1,Paneled->False]

System Models

sm=CreateSystemModel["sys1",{x'[t]- x[t]==u[t], y[t]==x[t]},t,{"u"∈"Modelica.Blocks.Interfaces.RealInput","y"∈"Modelica.Blocks.Interfaces.RealOutput"}]

sspec=<|"InputModel"->sm,"SamplingPeriod"->1|>;

cost=<|"Horizon"->3,"StateWeight"->{{5}},"InputWeight"->{{1}}|>;

cons=-2<=(Subscript[\[FormalX], 1])<=2∧-1<=(Subscript[\[FormalU], 1])<=1;

cd1=ModelPredictiveController[sspec,cost,cons]

cd1[{"ClosedLoopSystem",<|"ModelName"->"csys1"|>}]

Discrete input output model

MPC controller

cd=SystemsModelControllerData[];

cd["ClosedLoopSystem"]

List@@%

diom=cd["FeedbackGainsModel"]

diom["Path"];
%/.{n_, v_}:>{n,Head[v]->Union@Cases[v, e_[i_?IntegerQ]:>e[i],∞]}

diom=cd["QuasiFeedbackGainsModel"]

diom["Path"];
%/.{n_, v_}:>{n,Head[v]->Union@Cases[v, e_[i_?IntegerQ]:>e[i],∞]}

A general input output model

diom=DiscreteInputOutputModel[{u[0],u[1],u[2]},u]

OutputResponse[diom,{a,b,c,d,e}]

diom=DiscreteInputOutputModel[{(u[-1]+u[0])/2},u]

OutputResponse[diom, {a,b,c,d,e}]

diom=DiscreteInputOutputModel[{u[0]+y[-1]+y[-2],-y[0]+u[-2]},{{u,{-2,-1,1}}},{{y,{-1,-1,-1}}}]

or=OutputResponse[diom,Range[-20,20]];
ListStepPlot[or,DataRange->{0,40}]

Transfer function models

tfm=TransferFunctionModel[((z^(-1))/(0.1 (z^(-1))+1)),(z^(-1)),SamplingPeriod->1]

diom=DiscreteInputOutputModel[{ u[-1]-0.1 y[-1]},u,y ]

inp=Table[Sin[t],{t,0,2 π,0.1}];
GraphicsRow[Table[ListStepPlot[OutputResponse[sys,inp],
DataRange->{0,2π}],{sys,{tfm,diom}}]]

Properties

diom=DiscreteInputOutputModel[];

propsBasic={"Type","InputVariables","InputsCount","OutputVariables","OutputsCount","StatesCount","SamplingPeriod"};

Dataset[Thread[{propsBasic,diom[propsBasic]}],Alignment->Left]

propsCore={"Values","FirstInstant","FirstValue","Instants","LastInstant","LastValue","Path","PathComponent","PathComponents","PathFunction","PathLength"};

Dataset[Thread[{propsCore,diom[propsCore]}],Alignment->Left]

propsTime={"TimeValues","Times","FirstTime","LastTime","TimePath","TimeSeries","TemporalData"};

Dataset[Thread[{propsTime,diom[propsTime]}],Alignment->Left]

diom["Dataset"]

In closing...

Summary

• MPCs are computed based on a model, cost, and constraints

• MPC designs  are solved as parametric optimization problems

• The controller is in a closed form ready for further analysis and deployment

• The workflow is streamlined

Future work

• MPC speedup

• Control horizon

• Explicit tracking of nonconstant signals

• Full support of sampled-data systems in WL and discrete-time systems in WSM

Related resources

• ModelPredectiveController [ref]
• DiscreteInputOutputModel [ref]
• Wolfram R&D Livestreams [link]