Electricity prices.m

%% MPC

clear all
clear classes
clc

% User-defined parameters:
load('dwn_big.mat');
f0=0.1;
k=3000;
P.Hp=24;
P.Hu=P.Hp-1;
Hsim=100; % Simulation horizon >2
P.Wa=10;
P.Wx = P.Wx/1e6;

x0 = S.xmin + f0*(S.xmax-S.xmin);
epsilon=1e-3;

%svm data
Dd = DemandData(:,1); % My data
[yData,fData,sc_prms,Q,yscale]=svm_data(Dd,8760,k,200,P.Hp,1,1);
%load('mdl/svm_model1');   % load the svm model
load('eps_demand');   %load the epsion for the forecast
eps=eps./yscale;

%arima model
load('mdl/B2.mat');

% The closed-loop with the MPC controller:
CompTimes = zeros(Hsim-1,2); % 1st col: mpc_matrices time, 2nd column: computation time
mpc_ops = struct('use_soft_constraints',true,'use_beta',false, 'use_sparse', true,'use_forecast',true);
Pr_sim=4;
Xmpc=cell(Pr_sim,1);
Umpc=cell(Pr_sim,1);
%Xmpc=zeros(S.nx, Hsim);
%Umpc=zeros(S.nu, Hsim-1);

x=x0;
uMPC=[];

%Electricity cost
Eprice=cell(Pr_sim,1);
Eprice{1,1}=P.alpha2;
EconCost=cell(Pr_sim,1);
SmoothOpCost=cell(Pr_sim,1);
SafetyStorageCost=cell(Pr_sim,1);

Xmpc{1}(:,1)=x0;
err_percent=[0 0.001 0.005 0.01];
for i=1:Pr_sim
    Xmpc{i}(:,1)=x0;
    Eprice{i,1}=P.alpha2+err_percent(i)*randn(size(P.alpha2));
end

%% Compare all demands...
%load dwn_big
std_demand=zeros(S.nd,1);
range_demand=zeros(S.nd,1);
ele=zeros(S.nd,1);
for i=1:S.nd
    temp=(DemandData(:,i)./DemandData(:,1));
    std_demand(i,1)=std(temp);
    range_demand(i,1)=max(temp)-min(temp);
    ele(i,1)=DemandData(1,i)./DemandData(1,1);
end

%%
pastDataLength=k-5;

model='actual'; %selction of demand forecast method
dmd_err=zeros(S.nd*P.Hp,2);
%%
ctrlops=struct('normalise', false, 'mode', 'batch','u0', []);
for kk=1:Pr_sim
    P.alpha2=Eprice{kk,1};
    for j=k:k+Hsim-2
        if (~mod(j-k+1,10))
            disp(['Iteration : ' num2str(j-k+1)]);
        end
        switch model
            case 'svm'
                disp('svm');
                Q.Nt=j-201;
                Ypred=svmforecast(fData,yData,m,Q,1,sc_prms,yscale);
                Ypred=(Ypred./yscale)';
                D=vec(kron(ele', Ypred)'); % Choose D=Dreal to let the controller know the future..
                dmd_err(:,1)=vec(kron(ele', eps(:,1))'); % Choose D=Dreal to let the controller know the future...
                dmd_err(:,2)=vec(kron(ele', eps(:,2))'); % Choose D=Dreal to let the controller know the future...
            case 'arima'
                disp('arima');
                Dpast1 = DemandData(j-pastDataLength+1:j,1);
                Ypred = forecast(mdl.fit,P.Hp,'Y0',Dpast1);
                D=vec(kron(ele', Ypred)'); % Choose D=Dreal to let the controller know the future...
                dmd_err(:,1)=vec(kron(ele', eps(:,1))'); % Choose D=Dreal to let the controller know the future...
                dmd_err(:,2)=vec(kron(ele', eps(:,2))'); % Choose D=Dreal to let the controller know the future...
            otherwise
                disp('actual');
                D=vec(DemandData(j:j+P.Hp-1,:)'); % Future demand values - estimations/real?
        end
        Dreal=vec(DemandData(j:j+P.Hp-1,:)');
        tic;
        if (j==k)
            MPC = MPCmatrices(S, P, x, D, j, dmd_err,mpc_ops);
            controller = MPController(MPC, ctrlops);
        else
            %controller = controller.update(x, D, j, uMPC);
            controller = controller.update(x,D,j);
        end
        CompTimes(j-k+1,1)=toc;
        tic;
        [uMPC, Jstar, ystar, flag]= controller.control();
        %assert(flag==1 || flag==6,'CPLEX could not converge');
        if (flag==6)
            disp('WARNING : Non-optimal!');
            %disp(output);
        end
        assert(~isempty(uMPC), 'No control action could be computed');
        CompTimes(j-k+1,2)=toc;
        Useq=ystar(1:S.nu*(P.Hu+1),1);
        Xseq = MPC.Sx*x + MPC.Su*Useq + MPC.Sd*D;
        % Check whether the control action(s) satisfy the constraints and
        % whether the predicted sequence of states satisfies the state bounds.
        assert(max(Xseq-MPC.Xmax)<epsilon,'X>Xmax');
        assert(max(Xseq-MPC.Xmin)>epsilon,'X<Xmin');
        assert(max(Useq-MPC.Umax)<epsilon,'U>Umax');
        assert(max(Useq-MPC.Umin)>epsilon,'U<Umin');
        MPC = controller.getMPC();
        prInfeasibility=max(MPC.F ./ kron(  ones(1,size(MPC.F,2)),...
            abs(MPC.phi) )*ystar-sign(MPC.phi));
        disp(['Primal Infeas. : ' num2str(prInfeasibility)]);
        %assert(prInfeasibility<=1e-5,'Error: Inequality constraints not satisfied');
        % Check whether the state dynamics is followed:
        assert(norm(MPC.G*ystar-MPC.gamma)<=epsilon,'Error: Equality constraints not satisfied');
        dk = Dreal(1:S.nd);
        x = S.A*x + S.B*uMPC + S.Gd*dk;
        % Check whether the current state satisfied the constraints
        assert(sum(x> epsilon*ones(S.nx,1)+S.xmax)==0,'Critical Error: x>xmax');
        assert(sum(x< -epsilon*ones(S.nx,1)+S.xmin)==0,'Critical Error: x<xmin');
        assert(sum(uMPC<S.umin-epsilon*ones(S.nu,1))==0,'Critical Error: u<umin');
        assert(sum(uMPC>S.umax+epsilon*ones(S.nu,1))==0,'Critical Error: u>umax');
        Xmpc{kk}(:,j-k+2)=x;
        Umpc{kk}(:,j-k+1)=uMPC;
    end
    for i=1:Hsim-1
        EconCost{kk}(i,1)=P.Wa*(P.alpha1'+P.alpha2(k+i-1,:))*Umpc{kk}(:,i);
    end
    DeltaU = zeros(S.nu, Hsim-2);
    SmoothOpCost{kk}= zeros(1, Hsim-2);
    for i=1:Hsim-2
        DeltaU(:,i)=Umpc{kk}(:,i+1)-Umpc{kk}(:,i);
        SmoothOpCost{kk}(i) = DeltaU(:,i)'*P.Wu*DeltaU(:,i);
    end
    SafetyStorageCost{kk} = zeros(1, Hsim);
    for i=1:Hsim
        z=max(P.xs-Xmpc{kk}(:,i),zeros(S.nx,1));
        SafetyStorageCost{kk}(i)=z'*P.Wx*z;
    end
    SafetyStorageCost{kk} = SafetyStorageCost{kk}/1e7;
end
disp('All assertions were fulfilled!');
%% Plot the response:
% Plots of the response (repletion rate of the tanks) and the corresponding
% water demand.
for i=1:Pr_sim
    figure(1)
    plot(EconCost{i,1}(2:end));
    title('Electricity cost');
    hold all;
    figure(2)
    plot(SmoothOpCost{i,1}(2:end));
    title('Smooth operation cost');
    hold all;
    figure(3)
    plot(SafetyStorageCost{i,1}(2:end));
    title('SafetyStorageCost');
    hold all;
     average_cost(i,:)=[mean(EconCost{i,1}(2:end)) mean(SmoothOpCost{i,1}(2:end)) mean(SafetyStorageCost{i,1}(2:end))];
end 
figure(1)
legend('actual','0.1% error','0.5% error','1% error')
figure(2)
legend('actual','0.1% error','0.5% error','1% error')
figure(3)
legend('actual','0.1% error','0.5% error','1% error')

stat_to_plot=1;
for i=1:Pr_sim
    figure(4)
    title('State of the system')
    plot(Xmpc{i}(stat_to_plot,:))
    hold all;
end
legend('actual','0.1% error','0.5% error','1% error')

MPController.m

classdef MPController
    
    properties
        MPC_=MPCmatrices;
        uMPC_=[];
    end
    
    
    methods (Access=public)
        function obj = MPController(varargin)
            %controller = MPController(MPC);
            if isempty(varargin)
                return; % return empty object
            else
                obj=MPController();
            end
            obj.MPC_=varargin{1};
        end
        
        function [uMPC, Jstar, ystar, flag, output, time] = equ_ele_control(obj) %controller eliminating the equality constraints
            %g=null(full(obj.MPC_.G));
            %zo=obj.MPC_.G\obj.MPC_.gamma;
            G=full(obj.MPC_.G);
            n=size(G');
            [Q,R,P]=qr(G');
            R=R(1:n(2),1:n(2));
            Y1=Q(:,1:n(2));
            g=Q(:,n(2)+1:n(1));
            X_y1=R'\(P'*obj.MPC_.gamma);
            %H=full(obj.MPC_.H);
            zo=Y1*X_y1;
            assert(norm(G*zo-obj.MPC_.gamma)<10^-5,'not particular solution');
            f=obj.MPC_.f*g+zo'*obj.MPC_.H*g;
            H=sparse(g'*obj.MPC_.H*g);
            %size(H)
            Q=0.5*(H+H');
            %size(Q)
            %size([obj.MPC_.phi-obj.MPC_.F*zo;zo-obj.MPC_.ymin;obj.MPC_.ymax-zo])
            %size(g)
            %size(obj.MPC_.F*g)
            %size(sparse([sparse(obj.MPC_.F*g);-g;g]))
            l=sparse([sparse(obj.MPC_.F*g);-g;g]);
            tic;
            [ystar1, Jstar, flag, output]= ...
                cplexqcp(Q, f, l, ...
                [obj.MPC_.phi-obj.MPC_.F*zo;zo-obj.MPC_.ymin;obj.MPC_.ymax-zo], [], [], ...
                [], [], [],[], [],[]);
            time=toc;
            size(ystar1)
            size(l)
            S = obj.MPC_.getSystem();
            ystar=g*ystar1+zo;
            uMPC = ystar(1:S.nu,1);
        end
        
        function [uMPC, Jstar, ystar, flag] = control(obj)
            %[ystar, Jstar, flag, output]= ...
            %    cplexqcp(obj.MPC_.H, obj.MPC_.f, obj.MPC_.F, ...
            %    obj.MPC_.phi, obj.MPC_.G, obj.MPC_.gamma, ...
            %    [], [], [], obj.MPC_.ymin, obj.MPC_.ymax,[]);
            
            S = obj.MPC_.getSystem();
            tic;
            model.Q=2*obj.MPC_.H;
            model.obj=obj.MPC_.f;
            model.lb=obj.MPC_.ymin;
            model.ub=obj.MPC_.ymax;
            model.A=sparse([obj.MPC_.G;obj.MPC_.F]);
            model.rhs=full([obj.MPC_.gamma;obj.MPC_.phi]);
            no_eq=size(obj.MPC_.G,1);
            %model.sense=char([61*ones(no_eq,1);60*ones(size(obj.phi,1),1)]);
            model.sense=[repmat('=',no_eq,1);repmat('<',size(obj.MPC_.phi,1),1)];
            
            params.outputflag=0;
            result=gurobi(model,params);
            control.control_time=toc;
            if ~strcmp(result.status,'OPTIMAL')
                %control.flag=0;
                flag=6;
                warning('WARNING : Non-optimal!');
            else
                %control.ystar=result.x;
                %control.Jstar=result.objval;
                %control.flag=1;
                flag=1;
                Jstar=result.objval;
                ystar=result.x;
            end
            uMPC = ystar(1:S.nu,1);
        end
        
        function obj = update(obj,x, D, j)
            obj.MPC_ = obj.MPC_.update(x, D, j);
        end
        
        function mpc = getMPC(obj)
            mpc = obj.MPC_;
        end
        
    end
end