Vectorization

[Paul Draheim]

Hey at all!

I do have some problems about vectorization. I'm very new in programming in MatLab and Yalmip, so I hope you'll help me.

In the unit commitment example in the wiki, there's the little sentence: "remember, this can be vectorized, we use for-loop only to make the model easier to understand"

But how can I vectorize the code for the minimum up and down times?

for k = 2:Horizon
 for unit = 1:Nunits
  % indicator will be 1 only when switched on
  indicator = onoff(unit,k)-onoff(unit,k-1);
  range = k:min(Horizon,k+minup(unit)-1);
  % Constraints will be redundant unless indicator = 1
  Constraints = [Constraints, onoff(unit,range) >= indicator];
 end
end

Thanks!

[Johan Löfberg]

Is that a generic question, or do you mean this particular piece of code is extra tricky to vectorize.

The comment was very general in the sense that one should never forget that there might be simple vectorizations in this type of code which might both simplify the code, and reduce YALMIP overhead. The code above probably requires some clever use of repmat etc, i.e., not trivial for a MATLAB rookie, and won't improve readability, and only be important for performance if you have very long horizons or many units.

[Paul Draheim]

I do in general have my problems with the vectorization, but indeed I'd like to know how to vectorize this piece of code. I used it in a bit modified way for my own unit commitment. It will also have a horizon from a year, that means 8760 time steps and more. So I'm thankful for every vectorization... Could you please help me with that piece of code?

[Paul Draheim]

Dear Mr Löfberg,

do you have a solution for me?

[Johan Löfberg]

On closer inspections, that part is not easily vectorizable in the classical sense. If performance is essential, some index-generating approach would be possible

onoffindicies = [];
indicatorindicies1 = [];
indicatorindicies2 = [];

for k = 2:Horizon
 for unit = 1:Nunits

  range = k:min(Horizon,k+minup(unit)-1);
  v = ones(1,length(range));
  onoffindicies = [onoffindicies sub2ind([Nunits Horizon],unit*v,range)];
  indicatorindicies1 = [indicatorindicies1 sub2ind([Nunits Horizon],unit*v,k*v)];
  indicatorindicies2 = [indicatorindicies2 sub2ind([Nunits Horizon],unit*v,(k-1)*v)];
 end
end
Constraints = [Constraints, onoff(onoffindicies) >= onoff(indicatorindicies1)-onoff(indicatorindicies2)];

ouch.

BTW, I'm just testing some performance improvements which could be relevant for your model. Contact me per email if you want, could be a good benchmark.

[Paul Draheim]

Dear Mr Löfberg,

I sent you an Email. I hope that joh...@isy.liu.se is the right one!

Thanks in anticipation,

Paul

[Evgeny Schnittmann]

Dear Mr Löfberg,

I just started programming with YALMIP and have a similar problem. I have wrote a constraint to implement a maximum down-time for a unit in the following form:

for k = maxdown+1:Horizon
    pause1=sum(onoff(1,(k-maxdown):(k-1)));
    Constraints = [Constraints, implies(pause1==0, onoff(1,k)==1)];
end

Next I would like to vectorize the code, to make the simulation more efficient. I have tried to vectorize it myself for a couple of days now, but was not able to produce a satisfying solution. I would be very thankful for every kind of help.

Thanks in advance,

Evgeny

[Johan Löfberg]

I would just brute-force it, and then improve from there

You appear to have something like

horizon = 10;

maxdown = 3;

A = [ones(1,maxdown) zeros(1,horizon-maxdown)];

for i = 1:horizon-maxdown-1

    A = [A;circshift(A(end,:),1)]

end

pause1 = A*onoff(:);

implies(pause1 == 0, inoff(1,maxdown+1:Horizon)'==1)

and that can of course probably be improved on several levels (beyond fixing the bugs and indexing errors)

[Evgeny Schnittmann]

Thanks for the fast answer.

Unfortunately im getting this error when running the Code:

Error using  *  (line 624)
Inner matrix dimensions must agree.

Error in kevin_test (line 65)
pause1 = A*onoff(:);

Besides that I dont fully understand, what your code is trying to do.

[Johan Löfberg]

I simply looked at what indicies you are using when creating the sums

for k = maxdown+1:Horizon

  (k-maxdown):(k-1)
end

ans =

     1     2     3


ans =

     2     3     4


ans =

     3     4     5


ans =

     4     5     6


ans =

     5     6     7


ans =

     6     7     8


ans =

     7     8     9

and then I create a matrix which picks out those elements, and create all the sums using matrix-vector multiplication, although there was a typo and should be

A = [ones(1,maxdown) zeros(1,horizon-maxdown)];
for i = 1:horizon-maxdown-1

    A = [A;circshift(A(end,:),[0 1])]
end


A =

     1     1     1     0     0     0     0     0     0     0
     0     1     1     1     0     0     0     0     0     0
     0     0     1     1     1     0     0     0     0     0
     0     0     0     1     1     1     0     0     0     0
     0     0     0     0     1     1     1     0     0     0
     0     0     0     0     0     1     1     1     0     0
     0     0     0     0     0     0     1     1     1     0

it assumes onoff has length horizon, otherwise you will have to extract the relevant part accordingly

[Johan Löfberg]

or "simply"

A = toeplitz([ones(1,maxdown) zeros(1,horizon-maxdown)],[1 zeros(1,horizon-maxdown-1)])'