fmincon - penalty function

Hi,

I have the following problem:

I have a data set with half-hourly data for gas production and electricity demand. These are now to be brought together via a CHP (with a certain capacity). The produced gas is to be converted into electricity to cover the electricity demand as well as possible. So I am looking for a power P at each point in time.

To do this, I have calculated the deviation between demand and power in an objective function e and then summed it up.

q = fmincon(@(x) obj_func(x,P,demand,z,P_max,V_level,V_initial,V_size,V_proz,production,BHKW_mode,d),x0,[],[],[],[],lb,ub);
function e = obj_func(x,P,demand,z,P_max,V_level,V_initial,V_size,V_proz,production,BHKW_mode,d)

e = (demand - P).^2;
e = sum(e) + sum(d);

This objective function is then to be minimized with fmincon. This is done and at any time I get : P = demand (why should it not be so).

But now my constraint comes into play.

The gas comes first into a storage and if gas is converted into electricity, the storage content becomes smaller by this quantity.

V_level(1) = V_initial + production(1) - P(1);
for i = 2:z
    V_level(i) = V_level(i-1) + production(i) - P(i);
end
V_proz = 100*V_level/V_size;

And this storage may never be fuller than 100 % and never emptier than 0 %. I tried to realize this with a penalty function, which is added to e.

for i = 1:z
    if V_proz(i) > 100
        d(i) = 100000;
    elseif V_proz(i) < 0
        d(i) = 100000;
    else
        d(i) = 0;
    end
end

But no matter how big I choose the penalty for non-compliance with this condition: Nothing changes in the result, it is not considered by the optimizer at all....

Does anyone have any idea where my error lies or how I can solve it differently?

Thank you for your help, Matthias

Answer 1

[Torsten]

As far as I can see, you could formulate your problem as

min: sum_i (abs(demand(i)-P(i)))

s.c.

Pmin <= P(i) <= Pmax
0 <= V_level(i) <= V_size
V_level(i) = V_level(i-1) + production(i) - P(i)

If the difference to your former objective

min: sum_i (demand(i)-P(i))^2

is acceptable, this can be formulated as a linear optimization problem in the unknown P and V_level.

Use linprog to solve.

[Matthias K]

Thanks for your help Torsten!

The difference to my former objective would definitely be acceptable.

But it is still not clear to me, how I can implement those constraints:

Pmin <= P(i) <= Pmax
0 <= V_level(i) <= V_size
V_level(i) = V_level(i-1) + production(i) - P(i)

[Torsten]

Set the x-vector in linprog to

x = (P(1),P(2),...,P(n),V_level(1),...,V_level(n))

Then your constraints can be written as

Aeq*x = b  (V_level equality constraints)
A*x <= b   (other inequality constraints)

if you choose Aeq and A appropriately.

Just take a look at the examples provided in the documentation of linprog.

[Matthias K]

Ah I see, that's a smart idea, to let the optimizer determine V_level directly...

Thanks, I'll try that!

[Torsten]

Setting

eta_i = abs(demand(i)-P(i))

your problem formulation so far should read with the x-vector enlarged to

x=(P(1),...,P(n),V_level(1),...,V_level(n),eta(1),...,eta(n)): 

min: sum_i eta_i

s.c.

-eta_i <= demand(i) - P(i) <= eta_i
Pmin <= P(i) <= Pmax
0 <= V_level(i) <= V_size
V_level(i) = V_level(i-1) + production(i) - P(i)

[Matthias K]

Thanks again Torsten!

Can you explain why you introduce eta here?

I stayed with fmincon for the moment and I just got my first reasonable result. I just have to double check, if I have implemented this correctly.

V_level(i) = V_level(i-1) + production(i) - P(i)

I think Aeq and beq should be this, but I need to double check.

The first results look quite good. I don't want to show it here however, cause the plots contain unpublished research data...

[Torsten]

Can you explain why you introduce eta here?

If you want to stick to fmincon and use the quadratic objective, you don't need the eta formulation.

It's only necessary for the formulation as a linear optimization problem.

If the problem dimension gets larger and the problem remains linear, you should consider switching to linprog.

[Matthias K]

Just to wrap things up: I tried a lot of different approaches in the last days (nonlinear constraints, Intlinprog etc.) and learnt a lot about Matlab. In the end the solver patternsearch (global optimization toolbox) turnt out to be best. It even solved my problem, as I originally formulated it, even though it was written pretty dilettantely...

Just a matter of finding the right penalty costs.

[Torsten]

And which of the inputs to patternsearch did you fill ?

A, b, Aeq, beq, lb, ub, nonlcon ?

From the description, patternsearch should turn out to be quite safe, but extremely slow.

[Matthias K]

Just lb and ub.

I went back to just determine the values for P. so: x=(P(1),...,P(n))

Patternsearch didn't get along well with my many linear equality constraints, because then it doesn't find an initial point for its search which complies with the l.e. constraints.

Answer 2

[Matt J]

Problems

(1) The objective function code you've shown appears to depend on everything except the unknown variables x

(2) Assuming V_proz was supposed to be one of the unknowns, your penalty d is a discontinuous function of it.

[Matthias K]

Hi Matt,

thanks for your answer.

(1) Sorry, the variable x is used to determine the P. This is the full objective function. I wasn't really sure, what has to go into this function, so I just added every variable, which appeares in the code of the objective function.

function e = obj_func(x,P,demand,z,P_max,V_level,V_initial,V_size,V_proz,production,BHKW_mode,d)   

P_min = P_max * 0.25;

idx = find(x < P_min);
x(idx)=0;
P = x;

V_level(1) = V_initial + production(1) - P(1);

for i = 2:z
    V_level(i) = V_level(i-1) + production(i) - P(i);
end

V_proz = 100*V_level/V_size;
d = zeros(z,1);

for i = 1:z
    if V_proz(i) > 100
        d(i) = 100000;
    elseif V_proz(i) < 0
        d(i) = 100000;
    else
        d(i) = 0;
    end
end


e = (demand - P).^2;
e = sum(e) + sum(d);

end

(2) Could you elaborate on that?

[Matt J]

In that case,

(1) fmincon expects your objective to take a single input argument (a vector), not many. If you have multiple unkowns, they are expected to be the elements of the vector.

(2) The objective and constraint functions must be differentiable, and therefore continuous. Therefore discontinuous things like this are illegal,

if V_proz(i) > 100
    d(i) = 100000;
elseif V_proz(i) < 0
    d(i) = 100000;
else
    d(i) = 0;
end

I am not sure why you are using a penalty function instead of an explicit constraint. If V_proz are independent unknowns, why not just use fmincon's lb, ub inputs to bound them?

[Matthias K]

(1) Ah, I see, thank you. I'm still a novice, as you can probably tell.

(2) The reason, why I didn't use lb and ub is this:

I'm searching for the power at every time step, which is what the vector x contains. Since V_proz depends only indirectly on it, I didn't see how I could constrain V_proz via lb and ub.

Also if this works, I want to include another constraint. Can I rephrase the penalty function to make it work, or is that just the wrong approach?

[Matt J]

I didn't see how I could constrain V_proz via lb and ub.

Then use one of the other constraint inputs, like A,b, or nonlcon. The 'con' in fmincon stands for constrained. If you have no constraints, it makes little sense to be using fmincon - you may as well use fminunc.

I can't really tell from your code which variables are unknowns and which are fixed problem data, but it looks like V_proz is a linear function A*x of x, so linear inequality constraints of the form A*x<=1000, A*x>=0 should be all you need.

Incidentally, another discontinuous operation in your code is this,

idx = find(x < P_min);
x(idx)=0;

It should be removed. Instead, you should just impose the lower bound x>=P_max*0.25 with the lb input argument and increase all production(i) to be at least 0.25*P_max prior to the optimization,

production=max(production,0.25*P_max)

Here, I'm assuming P_max and production(i) are fixed, known problem data.

[Matthias K]

Thanks again Matt. Are you sure, that V_proz is a linear function of x?

Because V_proz is calculated from this. And P basically is x.

for i = 2:z
    V_level(i) = V_level(i-1) + production(i) - P(i);
end
V_proz = 100*V_level/V_size;

I guess, this is also a discontinous operation, which doesn't comply with fmincon, does it?

The lower limit for x is 0, which I gave as a constraint. The operation you talked about, was supposed to set all values between 0 and P_min to 0. But maybe I can solve that differently.

[Matt J]

I'm sure of nothing because again you haven't clarified which variables are unknown and what aren't. But the code for V_proz are all just additions and subtractions involving x=P, so it is quite linear with respect to x alone.

[Matthias K]

All variables apart from P, V_level and V_proz, which are calculated from x, are known problem data.

Here is my full code, if you want to look into it...

clear all
clc

 file_data = 'Daten_neu2.xlsx';

%% Read data

% Read data
opts = detectImportOptions(file_data,'Sheet','Results');
opts.VariableUnitsRange = 'A2';
data = readtable(file_data,opts,'Sheet','Daten');

% Assign variables

demand = data{:,2};     %Residual load
t = datenum(data{:,1});
t = t-t(1,1);           %time vektor
z = numel(t);
production = data{:,4}; %gas production


P_max = 300;           % Max capacity [kW]
P_min = P_max * 0.25;  % Min capacity [kW]

V_initial = 50;     % start gas storage level  [%]

V_size = 6;         % Gas storage size [h]
V_size = V_size * P_max/2;  %Gas storage size [kWh]

V_initial = V_initial/100*V_size;

P = zeros(size(t,1),1); % power
V_level = zeros(size(t,1),1);   %gas storage level [kWh]
V_proz = zeros(z,1);        %gas storage level [%]
d = zeros(z,1);
x0 = ones(z,1)*100; % Starting vector for optimization

%% Constraints

lb = zeros(z,1);
ub = ones(z,1);
ub = ub*P_max;

%% Solver
 
q = fmincon(@(x) obj_func(x,P,demand,z,P_max,V_level,V_initial,V_size,V_proz,production,d),x0,[],[],[],[],lb,ub);

%% calculation
P = q;

idx = find(P < P_min);
P(idx)=0;

V_level(1) = V_initial + production(1) - P(1);

for i = 2:z
    V_level(i) = V_level(i-1) + production(i) - P(i);
end

V_proz = 100*V_level/V_size;

for i = 1:z
    if V_proz(i) > 100
        d(i) = 100000;
    elseif V_proz(i) < 0
        d(i) = 100000;
    else
        d(i) = 0;
    end
end


Diff = (demand - P).^2;
SSQ = sum(Diff) + d;

%% Plot results
tiledlayout(2,1)

ax1 = nexttile;
plot(ax1,t,demand,'r');
hold on
plot(t,production,'Color','k');
hold on
plot(t,P,'Color','b');
hold off

ax2 = nexttile;
plot(ax2,t,V_proz,'r');


%% Objective function

function e = obj_func(x,P,demand,z,P_max,V_level,V_initial,V_size,V_proz,production,d)   

P_min = P_max * 0.25;

idx = find(x < P_min);
x(idx)=0;
P = x;


V_level(1) = V_initial + production(1) - P(1);

for i = 2:z
    V_level(i) = V_level(i-1) + production(i) - P(i);
end

V_proz = 100*V_level/V_size;
d = zeros(z,1);

for i = 1:z
    if V_proz(i) > 100
        d(i) = 100000;
    elseif V_proz(i) < 0
        d(i) = 100000;
    else
        d(i) = 0;
    end
end

e = (demand - P).^2;
e = sum(e) + sum(d);

end

[Matt J]

Here where you currently have,

idx = find(x < P_min);
x(idx)=0;

Are you sure that this isn't really supposed to be,

idx = find(x < P_min);
x(idx)=P_min;

Otherwise, V_proz will have jump discontinuities near x(i)=P_min.

[Matthias K]

yes, x(idx) shoudl be 0;

If x(idx) = P_min, I would have jump discontinuities near x(i)=0, wouldn't I?

But thanks a lot, for your answers, I've realized a lot of mistakes I made. I'll try to formulate my problem a little different

[Matt J]

If x(idx) = P_min, I would have jump discontinuities near x(i)=0, wouldn't I?

No, assuming P_min>0, it would be flat at x(i)=0, which you should see if you plot your obejctive as a function of any particular x(i).