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

6 次查看（过去 30 天）
jana nassereddine
回答： Emmanouil Tzorakoleftherakis 2023-6-13
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

请先登录，再回答此问题。

回答（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.

类别

在 Help Center 和 File Exchange 中查找有关 Controller Creation 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by