Subtractions in objective function

[ChrGr]

Hi Johan,

i was extending the unit commitment example but instead for power plants I modeled cooler plants. I already introduced a capacity charge and capacity revenue.

Other than the normal electricity charge ($/kWh) which is penalizing overall electricity consumption, this charge is penalizing high peak energy usage ($/kW) per months (1/30).

There are two capacity charges in my problem: 

• the capacity charge for the electricity consumed by the cooler plants (increasing the objective function because it's a cost): + 1/30*max(sum(CON(:,:),1))*1e-3*0.5*8;
  
• the capacity charge revenue by selling cold energy to the consumers (decreasing the objective function because it's a revenue): - 1/30*max(demand)*1e-3*20;

I attached the two files (main script and up-/downtime constraints-function) and put the main script here:

clearvars

%% Load data

% Energy price

cost_energy =[0.1050;0.1050;0.1050;0.1050;0.1050;0.1050;0.1050;0.1050;0.1050;0.1050;0.1050;0.1050;0.1050;0.1050;...

           0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;...

           0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1705;0.1050;0.1050];

% Relative load data            

demand_src = [0.15; 0.157; 0.164; 0.158; 0.1521; 0.1721; 0.1922; 0.1833; 0.1743; 0.1719; 0.1694; 0.1706; 0.1717; 0.1684; 0.1651; 0.2091; 0.253; 0.2877; 0.3224;...

             0.4388; 0.5553; 0.6178;0.6803; 0.7301; 0.7798; 0.8374; 0.8949; 0.9013; 0.9076; 0.9417; 0.9757; 0.9843; 0.9929; 0.9964; 1; 0.9618; 0.9236; 0.8892; ...

             0.8548; 0.8465; 0.8382;0.6507; 0.4632; 0.3134; 0.1636; 0.1609; 0.1582; 0.1582];

peak_load = 13650*1e3;                                       % Peak load to scale up the relative load data, in W

demand_src = demand_src.*peak_load;                              % Absolute AC load data 

% Plant Parameters

GENmax = 0.75.*[3500;3500;3500;3500;1800;1200;1200]*1e3;

GENmin = GENmax.*0.25/0.75;

eta = [7;7;7;7;6;5;5];                              % Efficiency (Coefficient of Performance)

minup   = [5;5;5;5;4;3;3];                          % Minimum uptime in time steps

mindown = [4;4;4;4;4;3;3];                          % Minimum downtime in time steps

% General parameters

Nunits = length(GENmax);                           % Number of power plants

Horizon = max([minup; mindown]);                    % 2.5 hrs

t_total = length(demand_src);                       % 1 day

% Extend data for moving horizon

demand_src = repmat(demand_src,ceil(t_total/Horizon),1);

cost_energy = repmat(cost_energy,ceil(t_total/Horizon),1); 

% Variables

Nhist = max([minup;mindown]);

demand = sdpvar(Horizon,1);

HistoryOnOff = sdpvar(Nunits,Nhist,'full');

onoff = sdpvar(Nunits,Horizon);

GEN = sdpvar(Nunits,Horizon);

CON = sdpvar(Nunits,Horizon);

ep = sdpvar(Horizon,1);

PreviusGEN = sdpvar(Nunits,1);

%% Constraints and Objective

constraints = [];

objective = 0;

for th = 1:Horizon

    % Constraints

constraints = [constraints, sum(GEN(:,th),1) >= demand(th)];

    for p=1:Nunits

        constraints = [constraints, onoff(p,th).*GENmin(p) <= GEN(p,th) <= onoff(p,th).*GENmax(p)]; %#ok<*AGROW>

        constraints = [constraints, CON(p,th) == GEN(p,th)*(eta(p)).^(-1)]; 

    end

    constraints = [constraints, min_updowntime_constraints([HistoryOnOff onoff],minup)];

    constraints = [constraints, min_updowntime_constraints(1-[HistoryOnOff onoff],mindown)];

    

    % Objective

objective = objective + sum(CON(:,th))*1e-3*0.5*ep(th);

objective = objective + 1/30*max(sum(CON(:,:),1))*1e-3*0.5*8;

objective = objective - 1/30*max(demand)*1e-3*20;

end

%% Optimizer

Parameters = {HistoryOnOff, demand, ep, PreviusGEN};

Outputs = {GEN, onoff};

% ops = sdpsettings('cachesolvers','1','solver','baron');

ops = sdpsettings('cachesolvers','1','debug','1','showprogress','1');

Controller = optimizer(constraints,objective,ops,Parameters,Outputs);

oldOnOff = repmat([0;0;0;0;0;0;0],1,Nhist);

oldGEN = repmat([0;0;0;0;0;0;0],1,Nhist);

figure

hold on; grid on;

for i=1:length(GENmax(:,1))

    legendstring2{i} = [num2str(GENmax(i).*1e-3/0.75), ' kW'];

end

legendstring2{end+1} = ['Demand'];

%% Moving Horizon

for t = 1:t_total    % 1 day

    disp(num2str(t))

    

    demand_out = demand_src(t:t+Horizon-1);

    energy_p = cost_energy(t:t+Horizon-1);

    

    [solution,problem] = Controller{{oldOnOff, demand_out, energy_p, oldGEN(:,end)}};

    optimalGEN = solution{1};

    optimalOnOff = solution{2};

    

    hold off

    stairs([-Nhist+1:0 1:Horizon],[oldGEN optimalGEN]');

    hold on

    stairs(1:Horizon,demand_out,'k');

    axis([-Nhist Horizon 0 max(demand_src)*1.1]);

    drawnow

    

    % Shift history

    oldGEN = [oldGEN(:,2:end) optimalGEN(:,1)];

    oldOnOff = [oldOnOff(:,2:end) optimalOnOff(:,1)];

    % Save variables

GEN_save(:,t) = optimalGEN(:,1);

end

legend(legendstring2,'Location','northwest')

%% Results

figure

hold on

bar(GEN_save','stacked')

for i=1:length(GEN_save(:,1))

    legendstring{i} = [num2str(GENmax(i).*1e-3/0.75), ' kW'];

end

legend(legendstring,'Location','northwest')

Because the capacity charge for selling cold energy is a revenue (for me as operator, but a charge for the cold consumer), it is decreasing the objective function (red line of code).

If I run this problem, it gets infeasible after one step:

+ Solver chosen : CPLEX-IBM

+ Processing objective function

+ Processing constraints

1

+ Calling CPLEX-IBM

2

Output argument "varargout{2}" (and maybe others) not assigned during call to "optimizer/subsref".

Error in tester_movinghorizon_cc (line 86)

    [solution,problem] = Controller{{oldOnOff, demand_out, energy_p, oldGEN(:,end)}};

If I change the minus to a plus (which is economically wrong), the problem solves fine but with a wrong objective function.

How can I address this problem? Unfortunately, I need the capacity revenue for my problem and can not scrap it. (Solver: CPLEX 12.6.3, Matlab version: R2016b, latest yalmip version)

Best regards,

Christoph

[chrgr]

I also tried different solvers. Gurobi, Mosek and Baron (with cplex and ipopt) also fail:

In an assignment  A(:) = B, the number of elements in A and B must be the same.

Error in optimizer/subsref (line 304)

            x(ismember(self.orginal_usedvariables,self.model.used_variables)) = output.Primal;

Error in tester_movinghorizon_cc_working (line 75)

[Johan Löfberg]

it crashes as you start sending data with NaNs to the optimizer object when t=2

[chrgr]

Thank you for your reply.

Yes, the NaN is somehow created in the first step by the 'max' operator and the minus sign.

What's confusing me is, that this only happens with the 'max' operator. The 'mean' operator for example works perfectly fine.

I now implemented a workaround with the 'mean' function. It gives the same results the solution range I'm observing (I'm recalculating the objective value afterwards with the final optimizer results).

Is there any way to work around the max function to get an proper objective function?

Best regards,

Christoph

[Johan Löfberg]

It looks extremely weird to me that your onoff variable is an sdpvar. Sounds like something that should be binary

When you subtract max in objective, you create nonconvexities which forces YALMIP to do MILP modelling of the max operator. However, as demand is a parameter, it makes no sense to complicate matters by adding this term to the objective. Further, as you have no lower bound on demand which is in the max operator, the MILP model will be very poor numerically, which can lead to all kinds of numerical issues

Also, as the objective is bilinear for undefined parameters, you have to explicitly tell YALMIP that the resulting model will be a MILP by specifying the solver (gurobi/cplex/mosek)

[chrgr]

The onoff variable as sdpvar was a copy and paste error. The optimizer solved it regardless of the variable definition.

I now explicitly defined a solver (cplex) and added bounds on 'demand' and it solved the problem immediately.

I know that it is pointless to calculated the max of demand inside the optimizer. Since it's a given parameter, I can compute the max outside and pass it to the solver. That'd yield the same result.

But this is just a simplification of my main problem where this 'max'-problem occurred for optimization variables where I don't know the value in advance. But putting bounds on those variables also solved the infeasibility also for my main problem.

Thank you very much for your help and the great tool yalmip!

[Johan Löfberg]

Well, you don't even need to compute it outside and pass it to the solver. It is constant, and thus is irrelevant