Find minimum consecutives in an optimvar

I try to "find" the minimum consecutives 0 or 1 in an integer optimvar.

Following a "normal logical" array, I would search for it as follows:

A = [0 0 0 1 1 0 0 1 0 1 1 1 0 0 0].';

% Find the position of the 1s

index = find(A);

% Calculate the difference between the zeros and find the indexes at wich a new block starts.

% A new Block starts at the first "index" and where the difference is bigger than 1

iStart = find([1; diff(index) >1]);

% The end of the block is always the index before the start of a new block

iEnd = [iStart(2:end)-1; size(index,1)];

% Calculate the length of the Blocks

onLength = iEnd-iStart+1;

offLength = ([index(iStart); length(A)+1] - [0; index(iEnd(1:end))] -1);

And I wanted to use this to "find" / define the minimum of a optimvar. But unfortunately this does not work, which is why I set the maximum like this:

% Define variables

prob = optimproblem;

n_max_off = 3;

N = 90;

ds_on = optimvar('ds_on', N, 'LowerBound',0, 'UpperBound', 1, 'Type','integer');

prob.Constraints.ds_max_off = optimconstr(N-n_max_off);

% The sum of the maximum +1 must be more the 1

% For more than the maximum, at least one element must be true.

for i = 1:N-n_max_off-1

prob.Constraints.ds_max_off(i) = sum(ds_on(i:i+n_max_off+1)) >= 1;

end

% Compressed code without loop

index = (0:n_max_off) + (1:N-n_max_off).';

prob.Constraints.ds_max_off(1:N-n_max_off) = sum(ds_on(index(1:N-n_max_off,:)),2) >= 1;

But to define the minimum I do not know how I can do it.

In other words:

I want to define a minimum and maximum length of zeros and ones which are my boundary conditions / constraints.

The maximum is not a problem, but to define the minimum length is my problem. So minimum and maximum length are my constraints.

My objective is someting different I left for readability.

[Julian]

I have found and already implemented a way to formulate the minimum number of a runtime linearly.

From the paper: "A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem"

DOI: 10.1109/TPWRS.2006.876672

I plan to post the solution in the near future.

Answer 1

[Julian]

The best solution is to formulate the minimum number of consecutives with a linear approach like in the paper: "A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem" DOI: 10.1109/TPWRS.2006.876672

At first define the variable:

% Define variables

prob = optimproblem;

n_min_off = 4;

n_max_off = 10;

n_min_on = 2;

n_max_on = 3;

N = 90;

I use two optimvars:

ds_on = optimvar('ds_on', N, 'LowerBound',0, 'UpperBound', 1, 'Type','integer');

ds_start = optimvar('ds_start', N, 'LowerBound',0, 'UpperBound', 1, 'Type','integer');

The maximum of the consecutives is like I implemented it in the question:

% Compressed code without loop

index = (0:n_max_off) + (1:N-n_max_off).';

prob.Constraints.ds_max_off = sum(ds_on(index(1:N-n_max_off,:)),2) >= 1;

% the same for maximum on

%...

I search for the start but in the end I musst not forget to give "ds_start" a little weight in my minimization function:

prob.Constraints.ds_start = diff(ds_on) <= ds_start(2:end);

Afterwards I say, after "ds_start" is 1, the following "ds_on"s musst be on too:

for i = 2:N-n_min_on+1:

prob.Constraints.min_off(i) = ds_start(i) * n_min_on - sum(ds_on(i:i+n_min_on-1)) <= 0;

end

% of faster:

index = (1:n_min_on) + (1:N-n_min_on).';

prob.Constraints.min_on = ds_start(2:N-n_min_on+1) * n_min_on - sum(ds_on(index),2) <=0;

And it's very siilar for "min_off":

for i = n_min_off+1:N

prob.Constraints.min_off(i) = ds_start(i) * n_min_on + sum(ds_on(i-n_min_on:i-1)) <= n_min_on;

end

% or faster:

index = (0:n_min_off-1) + (1:N-n_min_off).';

prob.Constraints.min_off = ds_start(1:N-n_min_off) * n_min_off + sum(ds_on(index),2) <= n_min_off;

With this calculation I'm able to define a minimum on and off time for a appliance in an optimization.

I am also able to make it multidemensional if I have more appliances of this kind.

With this workaround it is also possible to have different kinds of "on" values.

Answer 2

[Matt J]

It's going to be a nonlinear integer objective function, which means the solver it's going to choose is ga. You should save yourself the hassle and just implement the optimization with ga() directly.

[Julian]

Thank you for your quick answer.

But how should I define my minimum consraints with ga?

[Torsten]

And your constraint is that you want to limit the maximum length of subsequent zeros in your binary array ? Or is this your objective : min (maximum length of zeros) ?

[Julian]

I have a minimum and maximum length of zeros and ones.

The maximum is not a problem, but to define the minimum length is my problem. So minimum and maximum length are my constraints.

My objective is someting different I left for the readability.

[Torsten]

x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,intcon)

% Assuming A = x:
function [c,ceq] = nonlcon(x)
  A = x;
  
  n_max_off = ...;
  n_min_off = ...;
  % Find the position of the 1s
  index = find(A);
  % Calculate the difference between the zeros and find the indexes at wich a new block starts.
  % A new Block starts at the first "index" and where the difference is bigger than 1
  iStart = find([1; diff(index) >1]);
  % The end of the block is always the index before the start of a new block
  iEnd = [iStart(2:end)-1; size(index,1)];
  % Calculate the length of the Blocks
  onLength = iEnd-iStart+1;
  offLength = ([index(iStart); length(A)+1] - [0; index(iEnd(1:end))] -1);
  c(1) = max(offLength) - n_max_off;
  c(2) = -min(offLength) + n_min_off;
  ceq = [];
end

If you want to stick to problem-based formulation of your problem, use "fcn2optimexpr".

[c,~] = fcn2optimexpr(@nonlcon,x);
prob.Constraints.ds_max_off = c <= 0;

[Matt J]

You could also use this FEX download,

https://www.mathworks.com/matlabcentral/fileexchange/78008-tools-for-processing-consecutive-repetitions-in-vectors

>> A = [0 0 0 1 1 0 0 1 0 1 1 1 0 0 0];

>> [~,~,runlengths]=groupLims(groupTrue(logical(A)),1)

runlengths =

2 1 3

>> min(runlengths)

ans =

1

[Julian]

@Matt J I appreciate your help, but unfortunately the find function does not work with optimvar. Otherwise I could have used my code in the top of the question. But I think it great that you wrote this FEX.

@Torsten Do I understand it right that you would design it as a solver based problem? But I have some other constraints too, so is it ok to desing this constraints like this change it into a problem based formulation and add all the other stuff?

The second thing is, again the find function. When I try it the problem based way, my Input is an optimvar and I can't use the find function and I wasn't able to implement it otherwise.

[Torsten]

I gave you the problem-based command:

n_min_on = ...;
n_max_off = ...;
[c,~] = fcn2optimexpr(@nonlcon,x,n_min_on,n_max_off);
prob.Constraints.ds_min_on_max_off = c <= 0;

together with the function

function [c,ceq] = nonlcon(x,n_min_on,n_max_off)
  A = x;
  % Find the position of the 1s
  index = find(A);
  % Calculate the difference between the zeros and find the indexes at wich a new block starts.
  % A new Block starts at the first "index" and where the difference is bigger than 1
  iStart = find([1; diff(index) >1]);
  % The end of the block is always the index before the start of a new block
  iEnd = [iStart(2:end)-1; size(index,1)];
  % Calculate the length of the Blocks
  onLength = iEnd-iStart+1;
  offLength = ([index(iStart); length(A)+1] - [0; index(iEnd(1:end))] -1);
  c(1) = max(offLength) - n_max_off;
  c(2) = -min(onLength) + n_min_on;
  ceq = [];
end

x has to be defined before as a vector of integer optimvars bounded by 0 and 1.

[Bruno Luong]

As previously noted, using fcn2optimexpr you can convert any function to an expression. So in problem based optimization there is no reason for limiting to the built-in supported expression functions.

prob = optimproblem;
N = 10;
A = optimvar('A', N, 'LowerBound',0, 'UpperBound', 1, 'Type','integer');
minseqlgtexpr = fcn2optimexpr(@minseqlgt, A, 'OutputSize',[1,1],'Analysis','off','ReuseEvaluation',false);
prob.Constraints.mincst = minseqlgtexpr >= 1; % adjust the bound
maxseqlgtexpr = fcn2optimexpr(@maxseqlgt, A, 'OutputSize',[1,1],'Analysis','off','ReuseEvaluation',false);
prob.Constraints.maxcst = maxseqlgtexpr <= 6; % adjust the bound
% setup other stuffs and solve prob
    
function lgt = seqlgt(B)
lgt = diff(find([true; diff(B(:)~=0); true]));
end
function minlgt = minseqlgt(A)
minlgt = min(seqlgt(A));
end
function maxlgt = maxseqlgt(A)
maxlgt = max(seqlgt(A));
end

[Torsten]

Example:

n_max_on = 5;
n_min_off = 3;
N = 10;
prob = optimproblem;
ds_on = optimvar('ds_on', N, 'LowerBound',0, 'UpperBound', 1, 'Type','integer');
[c,~] = fcn2optimexpr(@nonlcon,ds_on,n_max_on,n_min_off);
prob.Constraints.ds_max_on_min_off = c <= 0;
prob.Objective = -sum(ds_on);
opts = optimoptions("ga",Display="none");
x0.ds_on = [1 1 1 1 1 0 0 0 1 1];
show(prob)

  OptimizationProblem : 

	Solve for:
       ds_on
	where:
       ds_on integer

	minimize :
       -ds_on(1) - ds_on(2) - ds_on(3) - ds_on(4) - ds_on(5) - ds_on(6) - ds_on(7) - ds_on(8) - ds_on(9) - ds_on(10)


	subject to ds_max_on_min_off:
       arg_LHS <= zeros(1, 2)

       where:

             [arg_LHS,~] = nonlcon(ds_on, 5, 3);

	variable bounds:
       0 <= ds_on(1)  <= 1
       0 <= ds_on(2)  <= 1
       0 <= ds_on(3)  <= 1
       0 <= ds_on(4)  <= 1
       0 <= ds_on(5)  <= 1
       0 <= ds_on(6)  <= 1
       0 <= ds_on(7)  <= 1
       0 <= ds_on(8)  <= 1
       0 <= ds_on(9)  <= 1
       0 <= ds_on(10) <= 1

[sol,fval,exitflag] = solve(prob,x0,Options=opts)

sol = struct with fields:
    ds_on: [10×1 double]

fval = -4

exitflag = 
    SolverConvergedSuccessfully

[c,~] = nonlcon(sol.ds_on,n_max_on,n_min_off)

c = 1×2

-1 0

sol.ds_on.'

ans = 1×10

0 0 0 1 1 1 1 0 0 0

function [c,ceq] = nonlcon(A,n_max_on,n_min_off)
  % Find the position of the 1s
  index = find(A);
  if isempty(index)
    c(1) = Inf;
    c(2) = Inf;
    ceq = [];
  else
    % Calculate the difference between the zeros and find the indexes at wich a new block starts.
    % A new Block starts at the first "index" and where the difference is bigger than 1
    iStart = find([1; diff(index) >1]);
    % The end of the block is always the index before the start of a new block
    iEnd = [iStart(2:end)-1; size(index,1)];
    % Calculate the length of the Blocks
    onLength = iEnd-iStart+1;
    offLength = ([index(iStart); length(A)+1] - [0; index(iEnd(1:end))] -1);
    c(1) = max(onLength) - n_max_on;
    c(2) = -min(offLength) + n_min_off;
    ceq = [];
  end
end

[Matt J]

@Julian But I have some other constraints too,

IMO, there is no real reason to be using problem-based optimization, unless there are complicated linear constraints that you haven't yet mentioned. If everything is nonlinear and you don't have multiple optimvar vectors, there is no benefit to continuing with the problem-based framework.

[Julian]

@Torsten and @Bruno Luong thank you for your explanation, now I got it. :)

I used the fcn2optimexpr wrong and had a bug in my code, I did not have here.

Thank you for your good explanation and patience, but unfortunately my computer keeps calculating and does not come to any result. Either I have a bug in the implementation or I have to set it up differently, which I will get to next week.

@Matt J I haven't mentiont it, because all of my other 4 implementet parts of my optimization are working fine. In addition, I want to build my optimization variable with respect to the number of individual elements, so that I can have several or no elements from one cathegory. And with the problem based framework I have a clear advantage.

[Julian]

I would have a second approach to defining the minima, as I am dissatisfied with the calculation time so far. I was thinking of the binomial formula, so that my values automatically become zero when the derivative is zero, i.e. when the current and subsequent states have not changed.

The formula for minimum on is:

My thought is, if my value changes, the derivative is not equal to 0. Subsequently, the n next elements must all be equal to 0 or 1, and sum to 0 or n_min_on.

The values, which do not interest me, I remove as follows:

If my value does not change, my derivative is equal to 0 => The following elements do not matter.

If my value changes, my derivative is +1 or -1 and to filter out the unwanted, I add - or +1.

As code, I wrote it as follows:

dif = [0; diff(ds_on)];

index = (1:n_min_on) + (1:N-n_min_on).';

prob.Constraints.ds_min_on = optimconstr(N-n_min_on);

prob.Constraints.ds_min_on(2:N-n_min_on+1) = dif(2:N-n_min_on+1).*(dif(2:N-n_min_on+1)+1).*(n_min_on-sum(2*(ds_on(index)-0.5),2)) == 0;

The formulation is a cubic problem (ds_on^3). I thought that Matlab might prefer a cubic equation and calculate it faster, but unfortunately I get the error:

Problems with integer variables and nonlinear equality constraints are not supported.

Do you think it is also possible to convert this equation using "fcn2optimexpr" or it doesn't make sense, because Matlab is faster with the other option?

One speed advantage I thought might be that "ReuseEvaluation" is true

Or it would possibly have an advantage, but rather with a direct "if" query in "fcn2optimexpr" and thus only one query and "linear" function?

[Torsten]

From MATLAB documentation:

x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,intcon) or x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,intcon,options) requires that the variables listed in intcon take integer values.

Note

When there are integer constraints, ga does not accept nonlinear equality constraints, only nonlinear inequality constraints.

[Julian]

I changed it to <= instead of == and it is calculating. I hope it will calculate faster, but I have still the question what kind of "problem" is better for Matlab to calculate?

[Torsten]

At least it's necessary to either remove integer variables or equality constraints from your code.