Valley-filling optimization problem

Hello, i'm trying to simulate a distributed scheduler for ev (electric vehicle) charging. the desired behavior would shave peaks and fill valley, in other words minimize the peak-to-average ratio of the electric load on a test grid. ( such as in this paper )

So, basically, what I'm doing is this:

• send to each user the aggregate load of the rest of the grid
• initialize the charge scheduling vector randomly, this will be the starting point x0 of the IPM optimization
• optimize: i want to minimize the peak-to-average ratio, or equivalently, the cost (quadratic fcn) of the supplied energy

The constraints are the following:

• The sum of the charges scheduled between now and the departure moment of the vehicle is equal to the energy requirement (first line of Aeq)
• Past scheduling cannot be changed (so it's equal to past load, ecs.ev_load ), scheduling post-departure must be zero ( I , the rest of Aeq matrix)
• lower bound of energy consumption is zero
• upper bound is a given parameter, ecs.max_load

x0=ecs.schedule;
I=eye(ecs.timeslots);
I(ecs.currenttime:ecs.ev.departure,:)=[];
Aeq=[zeros(1,ecs.currenttime-1) -1*ones(1,ecs.ev.departure-ecs.currenttime+1) zeros(1,ecs.timeslots-ecs.ev.departure);I];
beq=[-1*ecs.ev.requested_chg*ecs.nslotsperhour;ecs.ev_load;zeros(ecs.timeslots-ecs.ev.departure,1)];
A=[];
b=[];
nonlcon=[];
lb=zeros(ecs.timeslots,1);
ub=ecs.max_load*ones(ecs.timeslots,1);
options = optimoptions('fmincon','Algorithm','interior-point','MaxFunEvals',30000);

The objective function par at each ev sums

• x: the scheduled load of the vehicle (the optimization variable)
• other_loads: the future behavior of every non-ev load (unrealistic), the loads of the rest of ev schedules

other_loads=ecs.game_l;    %l_n (all schedules in the grid, mine included)
other_loads(:,ecs.id)=[];  %l_-n (removed my own schedule)
schedule = fmincon(@par,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);

function f=par(x)
  tot=sum([x other_loads],2);
  f=max(tot)/mean(tot);
end

Here's my problem: the optimization doesn't work. I keep tinkering with the code, but I always get (almost) the same result, where the valleys aren't filled at all, and the behaviour is almost constant. I verified by debugging that the content of the variables is correct, maybe I'm doing something wrong with the optimizer? Thanks in advance

[Matt J]

It would help if you showed us the output of

 whos Aeq beq x0 lb ub

so we could see the dimensions of everything.

[Riccardo]

>whos Aeq beq x0 lb ub
Name        Size              Bytes  Class     Attributes

Aeq       196x288            451584  double              
beq       196x1                1568  double              
lb        288x1                2304  double              
ub        288x1                2304  double              
x0        288x1                2304  double

Answer 1

[Matt J]

The function par() is not a twice differentiable function of x (because of the max() operation). It is therefore outside the scope of fmincon. Maybe try ga() instead?

Because you're not showing your actual code, there could be other problems, too. I'm assuming that the posted code is not the actual code, because you never pass other_loads to par() in any way. It should be returning an error message.

[Riccardo]

Sorry, i didn't want to clutter the post with much code. par() is a nested funcion, and fetches data from the parent function. By the way, here's the code: Here's the class ecs, representing a single house:

An array of houses is used to build a neighbourhood:

And finally, simulations is run by this script:

EDIT: code moved to attachment

[Matt J]

To reduce clutter, I took the liberty of moving the code from your previous comment to an attachment. However, I don't see the true version of par() or the code used to call fmincon anywhere in there. That was the important part.

[Riccardo]

Well, fmincon is

• called inside the function schedule=ipm_game(ecs),
• which is called by game(ecs) in the second 'if',
• called by reschedule(ecs),
• called by step_ecs
• called by step_nbhood
• called by the constructor of nbhood

Basically, ecs has got all the variables needed for the thing

Answer 2

[Matt J]

It also looks likely that the solutions to the problem are non-unique. The objective function only depends essentially on max(tot) and mean(tot). Given a solution tot with a certain max and mean, it is easy to find other tot with the same max and mean.

Your linear constraints Aeq*tot=beq may place some restriction on this, but it's not clear without saying Aeq,beq that it's enough.

[Riccardo]

Thank you for your answer. Would you recommend me to choose another objective function, in order to fill the valley?

[Matt J]

Now that I see that you have 196 equality constraints, I'm less concerned about this issue. However, differentiability might still be an issue, as I mentioned here.

[Riccardo]

This (badly written, sorry) quadratic cost function:

function f=game_objective_function(x) %(29) function to be minimized
              f=sum(cost(sum([x other_loads],2)));
          end
          function f=cost(x) % C_h=3/1600(E)^2
              f=(3/1600)*x.^2/ecs.nslotsperhour;
          end

should be differentiable, right? But it didn't bring me decent results last time. Now I'm running it again with FinDiffRelStep=0.1, and it's taking way longer, will post results asap

[Riccardo]

Nope. Still not what I want, but I see that (maybe) there is some microscopic response to the peaks of the inelastic demand.. Trying again with a steeper cost function now

[Riccardo]

Quite surprisingly, after removing the (3/1600) term from the cost function (minimum should stay the same), this is the behavior, as if it was 'overreacting':

Most probably the different ev's can't see each others' schedules

[Matt J]

This (badly written, sorry) quadratic cost function:

If the function is indeed quadratic then yes, it is differentiable, but if you're moving to a quadratic cost, it becomes more sensible to use quadprog instead of fmincon.

Now I'm running it again with FinDiffRelStep=0.1, and it's taking way longer, will post results asap

Regardless of whether you switch to quadprog, quadratic functions have simple derivatives, so it should be unnecessary to use finite difference derivative calculations. In particular, in fmincon, you should be supplying your own analytical gradient calculation, using the GradObj option, instead of tinkering with finite difference derivatives.

Answer 3

[Alan Weiss]

Without looking at your code, only having read your description, it seems to me that you are probably running into the problems described in the documentation about optimizing a simulation or ODE. Basically, if you use an Optimization Toolbox solver to optimize a simulation, then you need to take larger-than-default finite difference steps for gradient estimation.

The other thing to try is the patternsearch solver. I would use patternsearch before trying a genetic algorithm solver because patternsearch is much faster, more reliable, and easier to tune. But patternsearch is slower than an Optimization Toolbox solver when both have appropriate options set.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Riccardo]

Thanks, you're talking about the DiffMinChange parameter, right?

options = optimoptions('fmincon','Algorithm','interior-point','MaxFunEvals',30000,'DiffMinChange',0.1);

[Alan Weiss]

Yes, either DiffMinChange possibly combined with increased DiffMaxChange, or else FinDiffRelStep.