how to get the cost function result from model predictive controller?

Hello everyone,

so i have a model predictive controller programming code, and it works, but they don't give me the cost function as a resut when I run it, they only give me the evaluated performance of the tests as a result,

the code is:

% continuous-time state-space matrices

A=[1 0;0 1];

B=[0.0936 0.0936 0 0.0936;0 0.0752 0 0];

C=[1 0;0 1];

D=[0 0 0 0; 0 0 0 0];

Plant.InputName = {'MV1','MV2','MV3','MD'}; % set names of input variables

Plant.OutputName = {'MO1','MO2'}; % set names of output variables

Plant.StateName = {'S1','S2'}; % set names of state variables

Plant.InputUnit = {'Watt', 'Watt', 'Watt', 'Watt'}; % set units of input variables

Plant.OutputUnit = {'%', '%'}; % set units of output variables

Plant.StateUnit = {'%', '%'}; % set units of state variables

% create state space plant model

plant = ss(A,B,C,D);

% Assign Input and Output Signals to Different MPC Categories

plant = setmpcsignals(plant,MV=[1 2 3],MD=4)

Ts = 0.1; % sampling time

plant = c2d(plant,Ts); % convert to discrete time

% Display Basic Plant Properties and Plot Step Response

%damp(Plant)

% Plot the open-loop step response

%step(Plant)

% Create Controller

mpcobj=mpc(plant)

mpcobj.Ts = 0.1;

mpcobj.PredictionHorizon = 15;

mpcobj.ControlHorizon = 3;

% View and Modify Controller Properties

% Display a list of the controller properties and their current values.

get(mpcobj)

%mpcobj.DisturbanceVariables(1).Name= Power generated from PV panels;

% Override default scale factors using specified spans

%Uspan = [..., ..., ..., ...];

%Yspan = [..., ...];

% for i = 1:3

% mpcobj.MV(i).ScaleFactor = Uspan(iMV(i));

% end

% mpcobj.MD.ScaleFactor = Uspan(iMD);

% for i = 1:2

% mpcobj.OV(i).ScaleFactor = Yspan(i);

% end

%mpcobj.ManipulatedVariables(1).ScaleFactor = 10;

mpcobj.History = datevec(now);

MV = struct(Min=-1,Max=1);

%constraints

mpcobj.MV(1).MinECR = 0;

mpcobj.MV(1).MaxECR = 0;

mpcobj.MV(2).MinECR = 0;

mpcobj.MV(2).MaxECR = 0;

mpcobj.MV(3).MinECR = 0.2;

mpcobj.MV(3).MaxECR = 0.2;

mpcobj.MV(1).RateMinECR = 0.2;

mpcobj.MV(1).RateMaxECR = 0.2;

mpcobj.MV(2).RateMinECR = 5;

mpcobj.MV(2).RateMaxECR = 5;

mpcobj.MV(3).RateMinECR = 0.2;

mpcobj.MV(3).RateMaxECR = 0.2;

%mpcobj.MV(1).RateMax = ;

%mpcobj.MV(1).RateMin = ;

%mpcobj.MV(2).RateMax = ;

%mpcobj.MV(2).RateMin = ;

%mpcobj.MV(3).RateMax = ;

%mpcobj.MV(3).RateMin = ;

mpcobj.MV(1).Min = -2500;

mpcobj.MV(1).Max = 6000;

mpcobj.MV(2).Min = 0;

mpcobj.MV(2).Max = 6000;

mpcobj.MV(3).Min = 0;

mpcobj.MV(3).Max = 6000;

%power demand

%LOAD CURTAILMENT

%donnees: a= 0 pour 5-t-8 otherwise a =1 ; b= 0 pour 5-t-9 otherwise b=1; Np=...;

Pfdg=50;

Pfz= 22; PLi=267; Psh= 313; Pac=403; Ptv=86; Pwm=1.2; Pmicro=0.33; Pnot= 13; Pwater=8.5;

t=[0:1:24];

a=zeros(size(t)); a(t<5 | t>8)=1; % define constant by t

b=ones(size(t)); b(t>17 & t<21)=0;

Pdem=Pfdg+Pfz+PLi+(a+b)*(Psh+Pac+Ptv+Pnot)+a.*b*(Pwater+Pwm) + Pmicro;

plot(t,Pdem)

mpcobj.MV(3).Target=Pdem;

%mpcobj.MV.Name

%mpcobj.MV.Units

% define time-varying constraints and weights over the prediction horizon, to shifts at each time step

%mpcobj.MV.RateMin = [-2; -1.5; -1; -1; -1; -0.5];

%mpcobj.MV.RateMax = [2; 1.5; 1; 1; 1; 0.5];

%mpcobj.W.OutputVariables = [0.1 0; 0.2 0; 0.5 0; 1 0];

%weights

Weights.OV=5;% Above average priority

Weights.MV=[0 0 0.2 0]; % This value allows the MVs to move away from their targets temporarily to improve OV tracking.

Weights.MVRate=[0 0.1 0 0];

%Review Controller Design

review(mpcobj)

%custom set estimator

setEstimator(mpcobj,'custom');

xplant=[0;0];

xmpc = mpcstate(mpcobj);

SOCpbref=75; % (en pourcentage)

SOCliref=85; % (en pourcentage)

t1=25; % en seconds

Cmaxpb=22020; % As (Ampersecond)

Cmaxli=6000; % As (Ampersecond)

Ipb=[60 65 70 75 65 60 55 50 50 50 55 60 65 70 75 75 75 75 70 60 60 55 65 60 50]; % en Ampere

Ili=[70 75 80 85 75 70 65 60 60 60 65 70 75 80 85 85 85 85 80 70 70 65 75 70 60]; % en Ampere

Cpb= sum(Ipb)/t1;

Cli= sum(Ili)/t1;

SOCpb=(Cpb/Cmaxpb);

SOCli=(Cli/Cmaxli);

Nc=3;

for t = 0:Nc

y = (plant.C)*xplant; % plant equations: output

YY(t+1,:) = y;

xmpc.Plant = [SOCpb, SOCli]; % state estimation

u = mpcmove(mpcobj,xmpc,[],[SOCpbref, SOCliref]); % output is not needed

UU(t+1,:) = u;

u = [u;0];

xplant = plant.A*xplant + plant.B*u; % plant equations: next state

end

%Plot the simulation results

% figure

% subplot(2,1,1)

% plot(0:Ts:5,YY)

% title('y')

% subplot(2,1,2)

% plot(0:Ts:5,UU)

% title('u')

% Custom Cost Function

Optimization.ReplaceStandardCost = false

Optimization.CustomCostFcn = @myCostFunction;

% custom cost function

% e=2;

% CCli=1661.9;

% Cyclesli=3000;

% Cli=100;

% Vdc=48;

% effli=94;

% costdegre=10^(-9);

function J1 = myCostFunction(X,U,e,data,params);

for k=1:Np

Np = data.PredictionHorizon;

% Residential Loads Cost Function

U3 = U(1:p,data.MVIndex(3));

J1= e*sum(U1 -Pdem).^2

% Energy Storage System Cost Function

U2 = U(1:p,data.MVIndex(2));

J2= (CCli*U2*Ts*effli)/(2*Cyclesli*Cli) + costdegre*sum(U2.^2)

end

end

and how can uncomment the comments for better results;

thanks in advance

Answer 1

[Emmanouil Tzorakoleftherakis]

Please take a look at the doc page of mpcmove. The Info output containts a field called Cost. You can use it to visualize how the optimal cost changes across time steps.