Problem-based Optimization: Constraint sum(x==1) >=1 doesn't work.

Hello all,

i have a few questions:

1) I would like to show for a problem-based optimization problem that my variable x, a matrix, must contain the integer 1, 2 and 3 at least once. I have tried it like this:

A = [1 2 3;
     0 3 0]; 

[rowA, colA] = size(A);

prob = optimproblem("Description","Topologyoptimization");

x = optimvar("x",rowA,colA,"Type","integer","LowerBound",0,"UpperBound",3);


prob.Constraints.min1 = sum(x == 1) >= 1;
prob.Constraints.min2 = sum(x == 2) >= 1;
prob.Constraints.min3 = sum(x == 3) >= 1;

Unfortunately it does not work and this error message appears:

Error using sum

Invalid data type. First argument must be numeric or logical.

Error in test_problembased_topo (line 22)

prob.Constraints.mindF = sum(x == 1) >= 1;

What can I change to implement the constraints?

2) In addition, I want to create two boundary conditions which say that the integers 1 and 3 must always be located next to an integer 2 in the matrix (like it is in the matrix A). How can I implement this as a boundary condition?

[Matt J]

2) In addition, I want to create two boundary conditions which say that the integers 1 and 3 must always be located next to an integer 2 in the matrix (like it is in the matrix A). How can I implement this as a boundary condition?

Does that mean they will be horizontal neighbors, or could they also be vertical/diagonal neighbors?

[b999]

The could be horizontal or vertical neighbours, but not diagonal.

Answer 1

[Matt J]

I think it would be better to formulate it like this:

A = [1 2 3;
     0 3 0]; 

[rowA, colA] = size(A);

k=zeros(3); 
k([2 4 6 8])=1; %Neighbor-counting kernel

C=func2mat(@(z) conv2(z,k,'same'), zeros([rowA,colA]));%From the File Exchange:https://www.mathworks.com/matlabcentral/fileexchange/44669-convert-linear-functions-to-matrix-form

col=@(z) z(:);


x = optimvar("x",[rowA,colA,3],"Type","integer","LowerBound",0,"UpperBound",1);

prob = optimproblem("Description","Topologyoptimization");
prob.Constraints.noOverlap = sum(x,3) <= 1;
prob.Constraints.minContent = sum(sum(x,1),2)>=1;
prob.Constraints.Neighbors1 =  col(x(:,:,1)) - C*col(x(:,:,2)) <=0;
prob.Constraints.Neighbors3 =  col(x(:,:,3)) - C*col(x(:,:,2)) <=0;

...
    
xsol=solve(prob).x;
A=xsol(:,:,1) + 2*xsol(:,:,2) + 3*xsol(:,:,3);

[b999]

Thank you for helping me with the problem.

When it compiles, an error comes. Do I need to download the "func2mat"? I am quite new to Matlab

I have a few more questions:

1) why did you set the UpperBound to 1? I thought it has to be 3, beacause the integers are from 0 to 3

2) i want to implement some objective functions. In this case for example to reduce the numer of zeros in the MAtrix and maximize the 1s and 3s, without loosing the boundary condition. How is this possible ?

3) To me it looks like you have adjusted your code to the given matrix A. Finally I want to be able to enter a matrix of any size and it should work. Is it possible to write your code in a way that it adapts to the given matrix size?

[Matt J]

Do I need to download the "func2mat"?

Yes. It will give you the matrix C representing the convolution operation shown. The convolution conv2(z,k,'same') counts the number of neighbors around each point in z.

1) why did you set the UpperBound to 1? I thought it has to be 3, beacause the integers are from 0 to 3

Not anymore. In my formulation x(:,:,1) is a binary map showing the locations of the 1's. x(:,:,2) is a binary map of the 2's. x(:,:,3) is a binary map of the 3's.

2) i want to implement some objective functions. In this case for example to reduce the numer of zeros in the MAtrix and maximize the 1s and 3s, without loosing the boundary condition. How is this possible ?

The number of 1s is sum(col(x(:,:,1)). The number of 3's is sum(col(x(:,:,3)). You can build your objective from those expressions.

Is it possible to write your code in a way that it adapts to the given matrix size?

It already does, but be mindful that the optimzation solver will eventually choke if you make the problem too large.

[b999]

1) Is it possible to have a objective function which minimizes the 0s instead of maximizing the 1s or 3s?

2) Is it possible in genareal to have two or more objective functions which are weighted with factors? Like that maximizing the 1s is more important than maximizing the 3s?

2 ) I extended the code with a bigger matrix A, an objective function and a constraint

A = [1 2 2 3;
    0 2 0 0;
    0 3 0 0]; 
[rowA, colA] = size(A);

k=zeros(3); 

k([2 4 6 8])=1; %Neighbor-counting kernel

C=func2mat(@(z) conv2(z,k,'same'), zeros([rowA,colA]));%From the File Exchange:https://www.mathworks.com/matlabcentral/fileexchange/44669-convert-linear-functions-to-matrix-form

col=@(z)reshape(z,[],1);

x = optimvar("x",[rowA,colA,3],"Type","integer","LowerBound",0,"UpperBound",1);

prob = optimproblem("Description","Topologyoptimization");

prob.Constraints.noOverlap = sum(x,3) <= 1;
prob.Constraints.minContent = sum(sum(x,1),2)>=1;
prob.Constraints.Neighbors1 =  col(x(:,:,1)) - C*col(x(:,:,2)) <=0;
prob.Constraints.Neighbors3 =  col(x(:,:,3)) - C*col(x(:,:,2)) <=0;
prob.Constraints.Neighbors13 =  col(x(:,:,1)) - C*col(x(:,:,3)) >=1; % 1 shouldn't be next to 3

f_1 = -sum(col(x(:,:,1)));
f_2 = -sum(col(x(:,:,2)));

prob.Objective = (f_1 +f_2);

xsol=solve(prob).x;
A=xsol(:,:,1) + 2*xsol(:,:,2) + 3*xsol(:,:,3);

But now the code doesn't work. It says : Index in position 3 exceeds array bounds. Index must not exceed 1.

[Matt J]

1) Is it possible to have a objective function which minimizes the 0s instead of maximizing the 1s or 3s?

You can count the number of zeros with: numel(A)-sum(x(:)).

2) Is it possible in genareal to have two or more objective functions which are weighted with factors? Like that maximizing the 1s is more important than maximizing the 3s?

Sure. Sum the weighted combination to produce a single objective. Or see https://www.mathworks.com/discovery/multiobjective-optimization.html

2 ) I extended the code with a bigger matrix A, an objective function and a constraint

The new constraint makes the problem infeasible, so xsol=[]. I don't think the constraint expresses what you think it does.

[b999]

1) Ah thx, it works.

2) I implemented it in this way:

f_1 = -sum(col(x(:,:,1)));
f_2 = -sum(col(x(:,:,3)));

w_1 = 10;
w_2 = 8;

I dont know if it is the best solution in this case...

prob.Objective = (w_1*f_1 + w_2*f_2);

3) okay I see the problem. But right know i don't know how to implement this constarint. Do you know a solution for this right now? I thought I can do it like the constarints before.

I am happy that you are helping me with this problem. My target is to implement such a code for a matrix optimization with 9 integers( which should be arranged) and more constraints and objective functions. Maybe I need some help then. Is it okay, when i ask my questions in this chat?

[Matt J]

3) okay I see the problem. But right know i don't know how to implement this constarint. Do you know a solution for this right now? I thought I can do it like the constarints before.

The only way I can think of is,

prob.Constraints.Neighbors13 =  col(x(:,:,1))'*C.'*C*col(x(:,:,3)) <=0.9;
                    %equivalent to col(x(:,:,1))'*C.'*C*col(x(:,:,3)) = 0
                    %since LHS is integer

[Matt J]

Maybe I need some help then. Is it okay, when i ask my questions in this chat?

It would be better to Accept-click the answer given for the problem as you posted it (assuming the answer solves it). For more general problems you can post a new question.

[b999]

If I use

prob.Constraints.Neighbors13 =  col(x(:,:,1))'*C.'*C*col(x(:,:,3)) <=0.9;

the solver changes from intlinprog into ga. Why does it happen?

And when i am implemneting this constraint, the 1s and 3s are still next to each other.

I have to say I still don't get your approach with the 3-dimensional x and your function C and col. It is hard for me to understand how you implement the the Constarints and e.g. why you use the C and col. I think it is because i am quite new to Matlab.

[Torsten]

the solver changes from intlinprog into ga. Why does it happen?

Because the constraint is nonlinear. Intlinprog cannot cope with nonlinear constraints.

I think it is because i am quite new to Matlab.

Why does everybody write that it is due to MATLAB ? It is because you have no experience with optimization in general.

[Matt J]

I have to say I still don't get your approach with the 3-dimensional x and your function C and col. It is hard for me to understand how you implement the the Constarints and e.g. why you use the C and col.

Do you understand what the convolution operation does? From this example, you should see that convolution of z with k counts the number of non-zero horizontal and vertical neighbors for each pixel in z,

k=zeros(3); 
k([2 4 6 8])=1 %Neighbor-counting kernel

k = 3×3

0 1 0 1 0 1 0 1 0

z=   [0     1     0     1
      0     0     1     0
      0     1     0     1]

z = 3×4

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

     
Cz=conv2(z,k,'same')

Cz = 3×4

1 0 3 0 0 3 0 3 1 0 3 0

the solver changes from intlinprog into ga. Why does it happen?

Because Neighbors13 is a nonlinear constraint, we no longer have an integer linear program.

However, I think I have found a way to formulate the constraint linearly:

A = [1 2 2 3;
    0 2 0 0;
    0 3 0 0]; 
[rowA, colA] = size(A);

k=zeros(3); 

k([2 4 6 8])=1; %Neighbor-counting kernel

C=func2mat(@(z) conv2(z,k,'same'), zeros([rowA,colA]));%From the File Exchange:https://www.mathworks.com/matlabcentral/fileexchange/44669-convert-linear-functions-to-matrix-form

col=@(z)reshape(z,[],size(z,3));

x = optimvar("x",[rowA,colA,3],"Type","integer","LowerBound",0,"UpperBound",1);
y = optimvar("y",size(x),"Type","integer","LowerBound",0,"UpperBound",1); 


prob = optimproblem("Description","Topologyoptimization");

prob.Constraints.noOverlap = sum(x,3) <= 1;
prob.Constraints.minContent = sum(sum(x,1),2)>=1;

prob.Constraints.yUpper= col(y)<=C*col(x);    %this makes y a binary map of pixels with at least 1 neighbor
prob.Constraints.yLower= C*col(x)/8 <= col(y);

prob.Constraints.Neighbors1 =  col(x(:,:,1)) - C*col(x(:,:,2)) <=0;
prob.Constraints.Neighbors3 =  col(x(:,:,3)) - C*col(x(:,:,2)) <=0;
prob.Constraints.Neighbors13 =  x(:,:,1) +  y(:,:,3) <=1; % 1 shouldn't be next to 3

f_1 = -sum(col(x(:,:,1)));
f_2 = -sum(col(x(:,:,3)));
w_1 = 10;
w_2 = 8;

prob.Objective = (w_1*f_1 + w_2*f_2);

xsol=solve(prob).x;

Solving problem using intlinprog.
LP:                Optimal objective value is -92.500000.                                           

Cut Generation:    Applied 3 Gomory cuts, 35 clique cuts,                                           
                   12 cover cuts, and 4 zero-half cuts.                                             
                   Lower bound is -78.000000.                                                       

Heuristics:        Found 3 solutions using total rounding,                                          
                   and 1 solution using ZI round.                                                   
                   Upper bound is -78.000000.                                                       
                   Relative gap is 0.00%.                                                          


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05.

Asol=xsol(:,:,1) + 2*xsol(:,:,2) + 3*xsol(:,:,3)

Asol = 3×4

1 1 2 3 2 1 1 0 1 1 2 1

[b999]

@Torsten

Why does everybody write that it is due to MATLAB ? It is because you have no experience with optimization in general

Yes of course, you are right, I am not experienced in both areas, matlab and optimization.

[b999]

@Matt J thank you so much for your help.

Know I added some more integers, (1 to 8) and added some more constarints.

All in all the optimized matrix now looks strange: Asol =

30 30 30 30

30 30 30 30

10 30 30 30

24 10 30 30

Could you help me with this problems. That would be so nice. Thx!

A= [0 5 1 0;
    3 2 2 7;
    4 3 2 8;
    4 6 2 6]; 

[rowA, colA] = size(A);

k=zeros(3); 
k([2 4 6 8])=1; %Neighbor-counting kernel

C=func2mat(@(z) conv2(z,k,'same'), zeros([rowA,colA]));%From the File Exchange:https://www.mathworks.com/matlabcentral/fileexchange/44669-convert-linear-functions-to-matrix-form
col=@(z)reshape(z,[],size(z,3));

x = optimvar("x",[rowA,colA,8],"Type","integer","LowerBound",0,"UpperBound",1);
y = optimvar("y",size(x),"Type","integer","LowerBound",0,"UpperBound",1); 

prob = optimproblem("Description","Topologyoptimization");

prob.Constraints.noOverlap = sum(x,8) <= 1;
prob.Constraints.minContent = sum(sum(x,1),2)>=1;
prob.Constraints.yUpper= col(y)<=C*col(x);    %this makes y a binary map of pixels with at least 1 neighbor
prob.Constraints.yLower= C*col(x)/8 <= col(y);

prob.Constraints.One_nextto_Two =  col(x(:,:,1)) - C*col(x(:,:,2)) <=0; % 1 nextto 2
prob.Constraints.Three_nextto_Two =  col(x(:,:,3)) - C*col(x(:,:,2)) <=0; % 3 nextto 2
prob.Constraints.Six_nextto_Two =  col(x(:,:,6)) - C*col(x(:,:,2)) <=0; % 6 nextto 2
prob.Constraints.Seven_nextto_Two =  col(x(:,:,7)) - C*col(x(:,:,2)) <=0; % 7 nextto 2
prob.Constraints.Eight_nextto_Two =  col(x(:,:,8)) - C*col(x(:,:,2)) <=0; % 8 nextto 2
prob.Constraints.Five_nextto_One =  col(x(:,:,5)) - C*col(x(:,:,1)) <=0; % 5 nextto 1
prob.Constraints.Three_nextto_Four =  col(x(:,:,3)) - C*col(x(:,:,4)) <=0; % 3 nextto 4
prob.Constraints.Four_nextto_Four =  col(x(:,:,4)) - C*col(x(:,:,4)) <=0; % 4 nextto 4, there should be always minimum two 4s netxt to each other, is it correct to implement it like this?

prob.Constraints.One_not_nextto_Three =  x(:,:,1) +  y(:,:,3) <=1; % 1 shouldn't be next to 3
prob.Constraints.One_not_nextto_Three =  x(:,:,1) +  y(:,:,4) <=1; % 1 shouldn't be next to 4
prob.Constraints.One_not_nextto_Six =  x(:,:,1) +  y(:,:,6) <=1; % 1 shouldn't be next to 6
prob.Constraints.One_not_nextto_Seven =  x(:,:,1) +  y(:,:,7) <=1; % 1 shouldn't be next to 7
prob.Constraints.One_not_nextto_Eight =  x(:,:,1) +  y(:,:,8) <=1; % 1 shouldn't be next to 8

% How to implement, that the 1s should be always at the edge of the matrix ( so only at this cells A = [ 1 1 1 1; 1 0 0 1; 1 0 0 1; 1 1 1 1] for a 4x4 matrix
% How to implement that the 2s must have two neighbours, which are not 0, 4 or 5




%f_1 = -(numel(A)-sum(x(:))); % Maximize the 0s
f_1 = (numel(A)-sum(x(:))); % Minimize the 0s
%f_1 = -sum(col(x(:,:,1))); % Maximize the 1s
f_2 = 0;% -sum(col(x(:,:,3))); % Maximize the 3s

w_1 = 10;
w_2 = 8;

fun = w_1*f_1 + w_2*f_2;



prob.Objective = fun;

xsol=solve(prob).x;

Asol = xsol(:,:,1) + 2*xsol(:,:,2) + 3*xsol(:,:,3) + 4*xsol(:,:,4) + 5*xsol(:,:,5) + 6*xsol(:,:,6) + 7*xsol(:,:,7) + 8*xsol(:,:,8);

disp(int64(Asol));

[Matt J]

Could you help me with this problems. That would be so nice. Thx!

Yes, I can help, but a decision must be made first on whether your originally posted problem has been resolved by my answer, and if so, you should Accept-click it.

[b999]

Oh yeah, sorry. I'm not so long here on the forum. Is done :)

I wrote down a few questions in the previous post, maybe you know a clever solution for them. Meanwhile I try to understand the code a bit more detailed.

[Matt J]

A= [0 5 1 0;
    3 2 2 7;
    4 3 2 8;
    4 6 2 6]; 

[rowA, colA] = size(A);

k=zeros(3); 
k([2 4 6 8])=1; %Neighbor-counting kernel

C=func2mat(@(z) conv2(z,k,'same'), zeros([rowA,colA]));%From the File Exchange:https://www.mathworks.com/matlabcentral/fileexchange/44669-convert-linear-functions-to-matrix-form
col=@(z)reshape(z,[],size(z,3));

x = optimvar("x",[rowA,colA,8],"Type","integer","LowerBound",0,"UpperBound",1);
y = optimvar("y",size(x),"Type","integer","LowerBound",0,"UpperBound",1); 

prob = optimproblem("Description","Topologyoptimization");

prob.Constraints.noOverlap = sum(x,3) <= 1; %<----Matt J fixed
prob.Constraints.minContent = sum(sum(x,1),2)>=1;
prob.Constraints.yUpper= col(y)<=C*col(x);    %this makes y a binary map of pixels with at least 1 neighbor
prob.Constraints.yLower= C*col(x)/8 <= col(y);

prob.Constraints.One_nextto_Two =  col(x(:,:,1)) - C*col(x(:,:,2)) <=0; % 1 nextto 2
prob.Constraints.Three_nextto_Two =  col(x(:,:,3)) - C*col(x(:,:,2)) <=0; % 3 nextto 2
prob.Constraints.Six_nextto_Two =  col(x(:,:,6)) - C*col(x(:,:,2)) <=0; % 6 nextto 2
prob.Constraints.Seven_nextto_Two =  col(x(:,:,7)) - C*col(x(:,:,2)) <=0; % 7 nextto 2
prob.Constraints.Eight_nextto_Two =  col(x(:,:,8)) - C*col(x(:,:,2)) <=0; % 8 nextto 2
prob.Constraints.Five_nextto_One =  col(x(:,:,5)) - C*col(x(:,:,1)) <=0; % 5 nextto 1
prob.Constraints.Three_nextto_Four =  col(x(:,:,3)) - C*col(x(:,:,4)) <=0; % 3 nextto 4
prob.Constraints.Four_nextto_Four =  col(x(:,:,4)) - C*col(x(:,:,4)) <=0; % 4 nextto 4, there should be always minimum two 4s netxt to each other, is it correct to implement it like this?

prob.Constraints.One_not_nextto_Three =  x(:,:,1) +  y(:,:,3) <=1; % 1 shouldn't be next to 3
prob.Constraints.One_not_nextto_Three =  x(:,:,1) +  y(:,:,4) <=1; % 1 shouldn't be next to 4
prob.Constraints.One_not_nextto_Six =  x(:,:,1) +  y(:,:,6) <=1; % 1 shouldn't be next to 6
prob.Constraints.One_not_nextto_Seven =  x(:,:,1) +  y(:,:,7) <=1; % 1 shouldn't be next to 7
prob.Constraints.One_not_nextto_Eight =  x(:,:,1) +  y(:,:,8) <=1; % 1 shouldn't be next to 8

% How to implement, that the 1s should be always at the edge of the matrix ( so only at this cells A = [ 1 1 1 1; 1 0 0 1; 1 0 0 1; 1 1 1 1] for a 4x4 matrix
 x.UpperBound(2:end-1,2:end-1,1)=0;
 
% How to implement that the 2s must have two neighbours, which are not 0, 4 or 5
prob.Constraints.Neighbors_of_Two = 2*col(x(:,:,2)) <= C*col( sum(x(:,:,[1,3,6,7,8]),3));



%f_1 = -(numel(A)-sum(x(:))); % Maximize the 0s
f_1 = (numel(A)-sum(x(:))); % Minimize the 0s
%f_1 = -sum(col(x(:,:,1))); % Maximize the 1s
f_2 = 0;% -sum(col(x(:,:,3))); % Maximize the 3s

w_1 = 10;
w_2 = 8;

fun = w_1*f_1 + w_2*f_2;



prob.Objective = fun;

xsol=solve(prob).x;

Asol = xsol(:,:,1) + 2*xsol(:,:,2) + 3*xsol(:,:,3) + 4*xsol(:,:,4) + 5*xsol(:,:,5) + 6*xsol(:,:,6) + 7*xsol(:,:,7) + 8*xsol(:,:,8);

disp(int64(Asol));

[b999]

Thank you for helping me. I have a few more questions :)

1 ) I've added one more constraint, that the 2s always have to be next to another 2.

prob.Constraints.Three_nextto_Two =  col(x(:,:,2)) - C*col(x(:,:,2)) <=0; % 2 nextto 2

My target is, that all 2s should be connected, so like this ( so the 2s should be in a continous "way"):

A= [0 5 1 0;
    3 2 2 7;
    4 3 2 8;
    4 6 2 6]; 

I thought of adding something like this, so that the 2 ist part of the sum:

prob.Constraints.Neighbors_of_Two = 2*col(x(:,:,2)) <= C*col( sum(x(:,:,[1,2,3,6,7,8]),3));

But then ( with the objective function f_1 = -sum(col(x(:,:,1)));) the optimized matrix is:

5 1 1 1

1 2 2 2

2 7 8 3

2 6 4 4

Here you can see, that that not all 2s are connected continuously. Therefore the underlined 7 (7) for example has to be a 2. That would be neccessary for my problem. But the two constraints I mentioned does not fullfil this requirement yet.

2) Sometimes, e.g. if I use the objective function f_1 = (numel(A)-sum(x(:))); the matrix is

8 2 2 2

2 2 2 1

2 6 7 5

2 3 4 4

Here it would be good to have a constraint that avoid that there are 2s arranged in such a way

2 2 or 2 2 2 etc.

2 2 2 2 2

So one of the underlined 2s (2) shouldn't be there. The 2s can be considered as something like roads and in this case this arrangement doesn't make sense.

3) Is there a possible way to say, that for example a 1 has to be in a specific cell in the matrix. e.g. if you have a 3x3 matrix that there has to be a 1 at A(1,2)? I tried it like this but it doesn't work:

prob.Constraints.FixedValueConstraint = x(1, 2, :) == 1;

4) As a further objective function it would be nice, if you could maximize a "distance" between two intergers, e.g. maximize the Distance between the 1s and 4s. or that e.g. the 3s near to the 1s, but not neigbours...I don't know if it is possible, but that would be also very nice.

5) I Tried to implement, that there should be no 7 in the matrix, even if i definied 8 dimensions.

prob.Constraints.No_Seven =  sum(x(:,:,7)) <=0; % no 7s in the Solution

But it doesn't work. Is it only possible to implement something like that, if i reduce the dimension of x? It would be easier if I have always the dimension of 8 and exclude a specific dimension, because it depends on the application how much diemensions ( 3 to 8) i need.

6) When i am trying a 5x4 matrix it takes a lot of time to optimize it, like you said. Is there a way to speed up the whole optimization?

I apologize for all the questions. I try to solve them myself, but currently it is very difficult for me. Thank you very much again for your help. I appreciate it very much.

[Matt J]

The 2s can be considered as something like roads and in this case this arrangement doesn't make sense.

The constraint that you would need is that C*col(x(:,:,2)) is both (a) bounded from above by 2 and (b) contains only two 1's.

Is there a possible way to say, that for example a 1 has to be in a specific cell in the matrix. e.g. if you have a 3x3 matrix that there has to be a 1 at A(1,2)?

You can force this with the upper and lower bounds as in my last comment, or you could do

prob.Constraints.FixedValueConstraint = x(1, 2, 1) == 1;

As a further objective function it would be nice, if you could maximize a "distance" between two intergers, e.g. maximize the Distance between the 1s and 4s. or that e.g. the 3s near to the 1s, but not neigbours

You can define any objective or constraint function that you want using fcn2optimexpr. However, if the objective is nonlinear, as I think this would be, it will use ga() which is less reliable than intlinprog().

I Tried to implement, that there should be no 7 in the matrix, even if i definied 8 dimensions.

Again, you could use upper and lower bounds. Or,

prob.Constraints.No_Seven =  sum(col(x(:,:,7))) <=0

When i am trying a 5x4 matrix it takes a lot of time to optimize it, like you said. Is there a way to speed up the whole optimization?

You have 320 unknowns at this point. I guess that's not a small number....

You could try to provide an initial guess to intlinprog,

x = intlinprog(f,intcon,A,b,Aeq,beq,lb,ub,x0)

[b999]

Thx for your help!

The constraint that you would need is that C*col(x(:,:,2)) is both (a) bounded from above by 2 and (b) contains only two 1's.

I don't quite get what you mean by "bounded above" and "contains two 1's". The 2s should be togeteher and start next to a one, but they don't have to termine with a one. And should it be one constraint oder two? Do you have a eaxmple what you mean by your two expressions?

prob.Constraints.No_Seven = sum(col(x(:,:,7))) <=0

When i am using this i get:

Index in position 3 exceeds array bounds. Index must not exceed 1. Error in... Atopo = xsol(:,:,1) + 2*xsol(:,:,2) + 3*xsol(:,:,3) + 4*xsol(:,:,4) + 5*xsol(:,:,5) + 6*xsol(:,:,6) + 7*xsol(:,:,7) + 8*xsol(:,:,8);

It seems like, that it doesn't fit in this way. Do you mean by upper and lower bounds, bounds, that x.upper and x.lower should be zero for the whole matrix, like this?

x.UpperBound(1:end,1:end,6) = 0;
x.LowerBound(1:end,1:end,6) = 0;

This also doesn't work it has the same error.

[Matt J]

I don't quite get what you mean by "bounded above" and "contains two 1's". The 2s should be togeteher and start next to a one, but they don't have to termine with a one.

It means C*col(x(:,:,2))<=2 and C*col(x(:,:,2)) must contain only two 1s. We're not talking about A here.

Index in position 3 exceeds array bounds. Index must not exceed 1.

You've seen this before. It happens when the problem is infeasible and returns empty.

[b999]

You've seen this before. It happens when the problem is infeasible and returns empty.

Is it then possible to implemnt your approach? Because in my opinion these two expressions are mutually exclusive:

prob.Constraints.minContent = sum(sum(x,1),2)>=1; % here i define, that i want to have every integer in my solution

prob.Constraints.No_Seven =  sum(col(x(:,:,7))) <=0 % with this expression i say that i don't want the 7 in my solution

I think i couldn't use both expressions in one code. I thought of adding someting to the first expressin which expludes the 7.

It means C*col(x(:,:,2))<=2 and C*col(x(:,:,2)) must contain only two 1s. We're not talking about A here.

Okay, misunderstanding by me. But anyways i don't get your approuch of "C*col(x(:,:,2))<=2" and "C*col(x(:,:,2)) must contain only two 1s".

I would implement your thoughts like this (the second constraint doesn't work):

prob.Constraints.Continous_Two = C*col(x(:,:,2)) <=2;
prob.Constraints.Continous_Two_1 = sum(C*col(x(:,:,2))==1, 'all') == 2; % this one doesn't work because of the ==1, but i don't know how to implement it in anaother way

would mean, that the 2s only can have two neighbours (first constraint). And that two 2s have only one neighbour (second constraint). But why do i need these two constraint? If I want continous 2s they could have more than two neighbours, e.g.:

A= [0 5 1 0;

3 2 2 7;

4 3 2 8;

4 6 2 6];

and there are not necessarily two 2s with only one neighbour. Maybe don't understand these two constarints right? But even if they are right, currently i don't know how to implement them.

[Matt J]

Because in my opinion these two expressions are mutually exclusive:

You should indeed adjust the first expression,

prob.Constraints.minContent = sum(sum(x(:,:,[1:6,8]),1),2)>=1;

But anyways i don't get your approuch of "C*col(x(:,:,2))<=2" and "C*col(x(:,:,2)) must contain only two 1s".

You said the 2s are supposed to form a "road". If so, it means that every 2 along the road must have exactly two 2s as neighbors. The exception are the 2s at the very ends of the road, which will have only one neighbor. So the constraint might be formulated,

prob.Constraints.Continous_Two_1 = sum(C*col(x(:,:,2))  == 2*sum(  col(x(:,:,2)) )  - 2;

[b999]

prob.Constraints.minContent = sum(sum(x(:,:,[1:6,8]),1),2)>=1;

This works, thanks a lot!

But when i am using this expression

prob.Constraints.Continous_Two = C*col(x(:,:,2)) <=2;
prob.Constraints.Continous_Two_1 = sum(C*col(x(:,:,2)))  == 2*sum(col(x(:,:,2)))  - 2;

i still get this error: Index in position 3 exceeds array bounds. Index must not exceed 1.

[Matt J]

Well, a guaranteed way to avoid it is to provide a feasible initial guess, assuming you know one.

sol = solve(prob,x0)

[b999]

I already tried it like this

x0(:,:,1) = A == 1;
x0(:,:,2) = A == 2;
x0(:,:,3) = A == 3;
x0(:,:,4) = A == 4;
x0(:,:,5) = A == 5;
x0(:,:,6) = A == 6;
x0(:,:,7) = A == 7;
x0(:,:,8) = A == 8;
x0 = double(x0);
xsol=solve(prob,x0,"Options",options).x;

But it didn't work out, so I thought that i can't have a initial guess as an input. This is the error:

Error using optim.internal.problemdef.ProblemImpl/solveImpl

The value of 'x0' is invalid. Initial point must be a nonempty struct.

I know that A fullfills the constraints. And in my opnion I create a Matrix x0 which should be the size of xsol...obviously I did a mistake in defining x0.

[Matt J]

You need a struct that provides initial values for all of the optimvars.

x0=(A==reshape(1:8,1,1,[]));
y0=C*col(x0)>0;

sol0.x = x0;
sol0.y = y0;

xsol=solve(prob,sol0,"Options",options).x;

[b999]

In this case this error appears:

Error using optim.internal.problemdef.ProblemImpl/solveImpl

The value of 'x0' is invalid. Each field of initial point must contain double arrays.

[Matt J]

x0=double(A==reshape(1:8,1,1,[]));
y0=double(C*col(x0)>0);

[Torsten]

Most probably

y0 = double(C*col(x0)>0);

C*col(x0)>0 is of type "logical", not "double".

[b999]

x0=double(A==reshape(1:8,1,1,[]));

y0=double(C*col(x0)>0);

Now i get a new warning: Warning: x0 is infeasible. Ignoring x0.

But I am pretty sure that my A fullfills the constarints. Does this warning mean, that x0 does not have the right type or that A does not fullfills all constarints?

[Torsten]

It means that the initial guess vector does not satisfy the constraints.

[Matt J]

Use the infeasibility command to test each of your constraints and see which one is violated.

[b999]

I understand the command like this:

inf_test_ = infeasibility(prob.Constraints.Eight_nextto_Two,A(1,1));

or like this

inf_test_ = infeasibility(prob.Constraints.Eight_nextto_Two,A);

But it seems like the A does not fit as a pt like in the example by Matlab. (Error using optim.problemdef.OptimizationConstraint/getValue Second argument must be a nonempty struct.). In this case A is the input of the user.

Another point would be if it is possible to but in alle Constraints at once, so that you do not have to write so many infeasability tests? Maybe like this

inf_test_ = infeasibility(prob.Constraints,...);

[Torsten]

p.x = x0;
p.y = y0;
% + all other optimvars you have defined
inf_test = infeasibility(prob.Constraints.Eight_nextto_Two,p)

[Matt J]

But it seems like the A does not fit as a pt like in the example by Matlab

As before, infeasibility expects your struct sol0, not A.

Another point would be if it is possible to but in alle Constraints at once

It's not possible, but it is possible to loop over the Constraints, e.g.,

for con=string(fieldnames(prob.Constraints))'
    
    infTest.(con)= any(  infeasibility(prob.Constraints.(con),sol0)>0  , 'all');
    
end

[b999]

It works, perfect. Thank you very much @Torsten and @Matt J

Now i am litte bit confused, because when i am using

A= [0 1 5 0; 0 2 0 0;7 2 3 4;8 2 3 4; 0 6 0 0]

A = 5×4

0 1 5 0 0 2 0 0 7 2 3 4 8 2 3 4 0 6 0 0

The infeasability test says, that

prob.Constraints.Only_one_Five_nextto_One  = sum(x(:,:,5)) == sum(x(:,:,1));

is not fullfilt. But obviously it is fullfilled in A.

And the reason why i am doing this is to get continous 2s. But when i am implementing the two constraints @Matt J showed before

prob.Constraints.Continous_Two = C*col(x(:,:,2)) <=2;
prob.Constraints.Continous_Two_1 = sum(C*col(x(:,:,2)))  == 2*sum(col(x(:,:,2)))-2;

then the second one is also infeasable, but also in this case, like you can see it in A, it should be feasable.

[Torsten]

What does A have to do with the three-dimensional matrix x ?

[b999]

I am testing my input matrix A for feasability

x0= double((A==reshape(1:8,1,1,[]))); % A is the input matrix
y0= double(C*col(x0)>0);
sol0.x = x0;
sol0.y = y0;


for con=string(fieldnames(prob.Constraints))'
    infTest.(con)= any(  infeasibility(prob.Constraints.(con),sol0)>0  , 'all');
end

[Torsten]

A= [0 1 5 0; 0 2 0 0;7 2 3 4;8 2 3 4; 0 6 0 0];
x0 = A==reshape(1:8,1,1,[]);
sum(x0(:,:,5)) 

ans = 1×4

0 0 1 0

sum(x0(:,:,1))

ans = 1×4

0 1 0 0

[Matt J]

prob.Constraints.Continous_Two_1 doesn't work. It will need some more thinking.

[Matt J]

I think maybe the [5,1] constraint needs to be,

prob.Constraints.Only_one_Five_nextto_One  = C*col(x(:,:,5)) <= 4 - 3*col(x(:,:,1));

[b999]

I thougt of this solution regarding 5 and 1. Until now it also works i think.

prob.Constraints.Only_one_Five_nextto_One  = sum(sum(x(:,:,5))) == sum(sum(x(:,:,1))); % numbers 1 = numbers 5
prob.Constraints.Only_one_Five_nextto_One_2 = col(x(:,:,1)) <= C*col( sum(x(:,:,[5]),3)); % only one 5 next to 1

In addition to the variable y i added the variable yy, which should mean the neighbour of y (or the neighbour of the neighbour of x)

When i am implementing it like this:

prob.Constraints.yUpper= col(y)<=C*col(x);    %this makes y a binary map of pixels with at least 1 neighbor
prob.Constraints.yLower= C*col(x)/8 <= col(y);

prob.Constraints.yyUpper= col(yy)     <= C*col(y);  
prob.Constraints.yyLower= C*col(y)/8 <= col(yy);

prob.Constraints.Three_2_nextto_Six  =  x(:,:,3) -  yy(:,:,6) <=0; % 6 should be two cells nextto 3

it works for the optimization. but the feasability test is always 1 for yyupper and yyLower ( so both constraints are not feasable), no matter which A i have as an input.

[Matt J]

I don't know what K is.

[b999]

Sorry, K should be C...

Regarding yy:

I I implemented yy to have the possibility to say, that two fields next to 3 should be a 6. For example:

[ 0 1 0 0

0 2 0 0

0 6 2 3]

And therefor i thought it would be good to define the neighbours (yy) of the neighbours (y) of x

[ yy

yy y yy

yy y x y yy

.... ]

So like sthis, that i can define which integers should be two cells away from x.

And for this I thougt, that I can implement it like this

prob.Constraints.yyUpper= col(yy)     <= C*col(y);  
prob.Constraints.yyLower= C*col(y)/8 <= col(yy);

[Matt J]

Well, understanding why something is infeasible requires that we see both A and the code you used to generate sol0.

[b999]

This is my code. I still haven't found a solution for the continous 2s. In this case the infeasability of the Constarints which are containing yy is one.

clear all;
A= [5 1 0 0; 0 2 3 4;7 2 3 4;8 2 2 0; 0 2 6 0];
disp(A);

[rowA, colA] = size(A);

k=zeros(3); 
k([2 4 6 8])=1; %Neighbor-counting kernel

C=func2mat(@(z) conv2(z,k,'same'), zeros([rowA,colA]));%From the File Exchange:https://www.mathworks.com/matlabcentral/fileexchange/44669-convert-linear-functions-to-matrix-form
col = @(z)reshape(z,[],size(z,3));

x = optimvar("x",[rowA,colA,8],"Type","integer","LowerBound",0,"UpperBound",1);
y = optimvar("y",size(x),"Type","integer","LowerBound",0,"UpperBound",1); 

% yy as neighbour of neighbour of x
 yy = optimvar("yy",size(x),"Type","integer","LowerBound",0,"UpperBound",1);

prob = optimproblem("Description","Topologyoptimization","ObjectiveSense","minimize"); % Objective function should be minimized

prob.Constraints.yUpper= col(y)<=C*col(x);    %this makes y a binary map of pixels with at least 1 neighbor
prob.Constraints.yLower= C*col(x)/8 <= col(y);

prob.Constraints.yyUpper= col(yy)     <= C*col(y);  
prob.Constraints.yyLower= C*col(y)/8 <= col(yy);


prob.Constraints.noOverlap = sum(x,3) <= 1; 
prob.Constraints.minContent = sum(sum(x,1),2)>=1;

prob.Constraints.Only_one_Five_nextto_One  = sum(sum(x(:,:,5))) == sum(sum(x(:,:,1))); % For every 1 should be one 5
prob.Constraints.Only_one_Five_nextto_One_2 = col(x(:,:,1)) <= C*col( sum(x(:,:,[5]),3)); % One 1 should have only one 5 as neighbour

prob.Constraints.Min_two_Three = sum(x(:,:,3), 'all') <= 2; % Minimum 2 3s

x.UpperBound(2:end-1,2:end-1,1)=0; % 1s always at the edge of the matrix

prob.Constraints.One_nextto_Two =  col(x(:,:,1)) - C*col(x(:,:,2)) <=0; % 1 nextto 2
prob.Constraints.Three_nextto_Two =  col(x(:,:,3)) - C*col(x(:,:,2)) <=0; % 3 nextto 2
prob.Constraints.Two_nextto_Two =  col(x(:,:,2)) - C*col(x(:,:,2)) <=0; % 2 nextto 2 % new constraint
prob.Constraints.Six_nextto_Two =  col(x(:,:,6)) - C*col(x(:,:,2)) <=0; % 6 nextto 2
prob.Constraints.Seven_nextto_Two =  col(x(:,:,7)) - C*col(x(:,:,2)) <=0; % 7 nextto 2
prob.Constraints.Eight_nextto_Two =  col(x(:,:,8)) - C*col(x(:,:,2)) <=0; % 8 nextto 2
prob.Constraints.Five_nextto_One =  col(x(:,:,5)) - C*col(x(:,:,1)) <=0; % 5 nextto 1
prob.Constraints.Three_nextto_Four =  col(x(:,:,3)) - C*col(x(:,:,4)) <=0; % 3 nextto 4
prob.Constraints.Four_nextto_Four =  col(x(:,:,4)) - C*col(x(:,:,4)) <=0; % 4 nextto 4, there should be always minimum two 4s netxt to each other, is it correct to implement it like this?

prob.Constraints.One_not_nextto_Three =  x(:,:,1) +  y(:,:,3) <=1; % 1 shouldn't be next to 3
prob.Constraints.One_not_nextto_Four =  x(:,:,1) +  y(:,:,4) <=1; % 1 shouldn't be next to 4
prob.Constraints.One_not_nextto_Six =  x(:,:,1) +  y(:,:,6) <=1; % 1 shouldn't be next to 6
prob.Constraints.One_not_nextto_Seven =  x(:,:,1) +  y(:,:,7) <=1; % 1 shouldn't be next to 7
prob.Constraints.One_not_nextto_Eight =  x(:,:,1) +  y(:,:,8) <=1; % 1 shouldn't be next to 8
prob.Constraints.One_not_nextto_One =  x(:,:,1) +  y(:,:,1) <=1; % 1 shouldn't be next to 1


prob.Constraints.Neighbors_of_Two = 2*col(x(:,:,2)) <= C*col( sum(x(:,:,[1,2,3,6,7,8]),3)); % 2s must have two neighbours

prob.Constraints.Neighbors_of_Six = 2*col(x(:,:,6)) <= C*col( sum(x(:,:,[2]),3)); % 6s must have two neighbours which are 2s
prob.Constraints.Neighbors_of_Four = col(x(:,:,4)) <= C*col( sum(x(:,:,[3]),3));



prob.Constraints.Three_2_nextto_Six  =  x(:,:,6) -  yy(:,:,3) <=0; % 6 should be two cells nextto 3


% Continous 2s --> does not work
% prob.Constraints.Continous_Two = C*col(x(:,:,2)) <=2;
% prob.Constraints.Continous_Two_1 = sum(C*col(x(:,:,2)))  == 2*sum(col(x(:,:,2)))-2;

%% Objective

f_1 = -sum(col(x(:,:,1))); % Maximize the 1s
f_2 = -sum(col(x(:,:,3))); % Maximize the 3s

w_1 = 1;
w_2 = 1;

fun = w_1*f_1 + w_2*f_2;

prob.Objective = fun;

options = optimoptions('intlinprog','Display','iter');


x0= double((A==reshape(1:8,1,1,[])));
y0= double(C*col(x0)>0);
yy0 = double(C*col(x0>0));
sol0.x = x0;
sol0.y = y0;
sol0.yy = yy0;


for con=string(fieldnames(prob.Constraints))'
    
    infTest.(con)= any(  infeasibility(prob.Constraints.(con),sol0)>0  , 'all');
    
end

xsol=solve(prob,sol0,"Options",options).x;


Asol = xsol(:,:,1) + 2*xsol(:,:,2) + 3*xsol(:,:,3) + 4*xsol(:,:,4) + 5*xsol(:,:,5) + 6*xsol(:,:,6) + 7*xsol(:,:,7) + 8*xsol(:,:,8);

disp(int64(Asol));

[Matt J]

yy0 = double(C*col(y0)>0);

[b999]

When i am using this i get

Error using *

Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform elementwise multiplication, use '.*'.

[Matt J]

Check that the number of columns in the first matrix matches the number of rows in the second matrix.

Indeed, you should check that.

[b999]

So C has a dimension of 20x20

x0: 5x4x4x8

y0: 20x8 and col(y0): 160x1

So C* col(y0) would be 20x20 * 160x1, which doesn't fit.

So I think it has to be

yy0 = double(C*(y0)>0);

--> Edit: seems to work.

Then I can concentrate now on the continuity of the 2s

[Matt J]

Edit: seems to work.

Yep. That does seem right.

[b999]

If got another question regarding the implementation of if-conditions.

For example I want number 1s = number of 5s, except the number of 1s is bigger than two.

I tried to implement it like that:

if sum(sum(x(:,:,1))) >= 3
    prob.Constraints.Only_one_Five_nextto_One  = sum(sum(x(:,:,5))) == 2*sum(sum(x(:,:,1)));
end

But the error is: Conversion to logical from optim.problemdef.OptimizationInequality is not possible.

Is it in general not possible to implement such if conditions in this kind of optimization?

I thought maybe it would be an opportunity somehow to implement among other things the continous 2s.

[Torsten]

Run the optimization two times:

Once with the constraint sum(sum(x(:,:,1))) >= 3 and with the above constraint and once with the constraint sum(sum(x(:,:,1))) < 3 and without the above constraint and see which simulation yields the better result.

[Matt J]

Regarding the continuity of the 2's, I believe you would want the following.

xc=col(x(:,:,2));

xc.'*C*xc == 2*sum(xc)-2;

Unfortunately, it is a quadratic constraint, so the solver will resort to ga. Moreoever, ga does not allow both integer and equality constraints, so you would have to relax it to,

 1.5 <= 2*sum(xc)-xc.'*C*xc <= 2.5;

Hopefully, a linear alternative can be found.

[b999]

@Torsten For one if-condition it might work. But if I have e.g. 3 if-conditions, which are similarly built, I would have to run eight different versions of the code. I think that wouldn't be a perfect solution.

@Matt J thanks for helping. I will check your solution. But i think a ga algorithm wouldn't be the best solutuon. Especially for a comparison of the optimizations.

I have thought about the boundary condition that the 2s should be arranged continuously. I think for my case it would be also okay, that, in the case that there are several 1s (or of course only one 1), the 2s must lead all to a 1. But not all 2s must be connected with all 2s. There must be no 2s without connection to a 1.

What i mean: this would be also okay regarding the 2s:

A_okay = [1 2 2 8;
     5 3 3 6;
     0 4 4 2;
     0 0 7 2;
     0 0 5 1]

A_okay = 5×4

1 2 2 8 5 3 3 6 0 4 4 2 0 0 7 2 0 0 5 1

But e.g. this wouldn't be okay:

A_notokay=   [ 1   2   2   8;
               5   3   3   6;
               3   4   4   2;
               2   3   3   2;
               2   7   5   1]

A_notokay = 5×4

1 2 2 8 5 3 3 6 3 4 4 2 2 3 3 2 2 7 5 1

The left 2s are not connectetd to a 1.

I am tyring to solve this. Maybe it is a llite bit easier. The 2s could be of course three or more in a row. The two matrices are only examples. @Matt J do you have maybe a first approach for such a condition?

I've got another question regarding saying that a 4 should be in the center of the matrix. I implemented this expression that the 4s should not be at the edge of the matrix. This seems to work

x.UpperBound([1:end],[1,end],4) = 0;   % No 4 at the edge of the matrix
x.UpperBound([1,end],[1:end],4) = 0;   % No 4 at the edge of the matrix

But the best solution would be if I could have a constraint wich says, that the 4 should be as far in the center of the optimized matrix as possible. And not only, that the 4s shouldn't be at the edge. It would be best if the condition does not contain any restriction of rows and columns, since you do not know how many terminals you will have in the end.

[b999]

Besides the two questions one post before, if got another question regarding of the understanding of this expressions:

prob.Constraints.yUpper= col(y)<=C*col(x);    %this makes y a binary map of pixels with at least 1 neighbor
prob.Constraints.yLower= C*col(x)/8 <= col(y);

The first constarint defines y y-matrix, which has at least one neighbor to every x, is it correct?

I didn't get the second constraint. It defines the minimum of y, but why do I divide the vector of the neighbours of x by 8? Why do I have do devide the vector of the neighbours and why should I divide it by 8 and not another number?

I would be very happy if you could help me with this question. And of course also with the question asked before. Thank you very much.

[Matt J]

yUpper ensures that when x(i,j,k) has no neighbors y(i,j,k) will be 0.

yLower ensures that when x(i,j,k) has a neighbor y(i,j,k) will equal 1.

The factor of 1/8 was not a unique choice. Any number ensuring that the left hand side will always be <=1 would suffice.

[Matt J]

Here's an idea for continuous 2s:

z=optimvar('z',[rowA,colA],'type','integer','Lower',0, 'Upper',2)

Constraints.conTwos_1=  z<=2*(x(:,:,2));

Constraints.conTwos_2=  z<= C*col(x(:,:,2));
Constraints.conTwos_3=  z>= C*col(x(:,:,2))-4*(1-x(:,:,2));

Constraints.conTwos_4=  sum(z(:))==2*sum(col(x(:,:,2)))-2;

The first constraint implies that z can only be non-zero where x(:,:,2) is nonzero. The next two constraints ensure that everwhere x(:,:,2)=1 will have z=C*col(x(:,:,2)). In other words, z indicates how any neighboring 2s each of the 2s have.

The fourth constraint, together with the first, ensures that all 2s have exactly two neighboring 2s, except for those at the tips of the strand, which will have only one neighboring 2.

[b999]

Thanks for your approach.

I think the dimensions doesn't fit everywhere ( e.g. dim(z) ~= dim(C*x)). I think you have to converte the z with col ( i don't know if the col is always at the right place)

z=optimvar('z',[rowA,colA],'type','integer','Lower',0, 'Upper',2);
Constraints.conTwos_1 =  z <= 2*(x(:,:,2));
Constraints.conTwos_2 =  col(z) <= C*col(x(:,:,2));
Constraints.conTwos_3 =  col(z) >= C*col(x(:,:,2)) - 4*col((1-x(:,:,2)));
Constraints.conTwos_4 =  sum(z(:))==2*sum(col(x(:,:,2)))-2;

But when i am using it like this, the solution is still e.g.:

A=[  6   2   1   5;
    3   2   2   8;
    4   7   3   4;
    4   3   3   4;
    3   2   2   3]

A = 5×4

6 2 1 5 3 2 2 8 4 7 3 4 4 3 3 4 3 2 2 3

Wouldn't you also have to have a condition that it starts at a 1? So something with x(:,:,1)?

I am referring to this:

I have thought about the boundary condition that the 2s should be arranged continuously. I think for my case it would be also okay, that, in the case that there are several 1s (or of course only one 1), the 2s must lead all to a 1. But not all 2s must be connected with all 2s. There must be no 2s without connection to a 1.

I think you are referring to the one which connects all 2s. This would be also possible, but maybe not the preferred solution for my problem. But maybe one can arrange this solution to a solution which connects the 2s to different 1s, as described before. I hope it is not too confusing.

[Matt J]

I think you are referring to the one which connects all 2s.

Yes. I thought that's what we were calling "continuous twos". The pattern of multiple strings of twos and ones would be a fair bit harder, I think.

But when i am using it like this, the solution is still e.g.:

Odd, I wasn't getting anything that bad, but in the (deliberately simplified) example below, I've improved things (and gotten rid of col():

[rowA, colA, maxA] = deal(8,6,5);

k=zeros(3); 
k([2 4 6 8])=1; %Neighbor-counting kernel
C=func2mat(@(z) conv2(z,k,'same'), zeros([rowA,colA]));


x = optimvar("x",[rowA*colA,maxA],"Type","integer","LowerBound",0,"UpperBound",1);



Constraints.noOverlap = sum(x,2) == 1;
Constraints.minContent = sum(x)>=1;
Constraints.twosUpper = sum(x(:,2))<=rowA+colA; %just to keep the 2s sparse

z=optimvar('z',rowA*colA,'Type','integer','Lower',0, 'Upper',2);
Constraints.conTwos_1 =  z <= 2*(x(:,2));
Constraints.conTwos_2 =  z <= C*x(:,2);
Constraints.conTwos_3 =  z >= C*x(:,2) - 4*(1-x(:,2));
Constraints.conTwos_4 =  sum(z)==2*sum(x(:,2))-2;
Constraints.conTwos_5 =  noBlocks(x,rowA,colA);


prob = optimproblem("Constraints",Constraints,'ObjectiveSense','max');
prob.Objective=sum(x(:,2));

sol=solve(prob);

Solving problem using intlinprog.
LP:                Optimal objective value is 14.000000.                                            

Cut Generation:    Applied 1 cover cut, 6 zero-half cuts,                                           
                   and 1 flow cover cut.                                                            
                   Upper bound is 14.000000.                                                        

Heuristics:        Found 3 solutions using diving.                                                  
                   Lower bound is 14.000000.                                                        
                   Relative gap is 0.00%.                                                          


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05.

A=round(reshape(sol.x*(1:maxA)', rowA,colA))

A = 8×6

4 4 5 1 3 3 5 2 2 2 5 3 4 4 5 2 4 5 1 5 2 2 4 4 4 2 2 5 4 5 5 2 4 5 4 5 1 2 5 2 4 3 1 2 2 2 4 1

A==2

ans = 8×6 logical array
   0   0   0   0   0   0
   0   1   1   1   0   0
   0   0   0   1   0   0
   0   0   1   1   0   0
   0   1   1   0   0   0
   0   1   0   0   0   0
   0   1   0   1   0   0
   0   1   1   1   0   0

function  con=noBlocks(x,rowA,colA)

    k0=zeros(3); k0([1,2,4])=1; 
    
    for i=4:-1:1
    
      k=rot90(k0,i);
      K=func2mat(@(z) conv2(z,k,'same'), zeros([rowA,colA]));
    
      con{i}= K*x(:,2)<=9*(1-x(:,2))+2*x(:,2);
    
    
    end
    con=[con{:}];

end

[Matt J]

The pattern of multiple strings of twos and ones would be a fair bit harder, I think.

Or maybe not? See below:

[rowA, colA, maxA] = deal(8,6,5);

k=zeros(3); 
k([2 4 6 8])=1; %Neighbor-counting kernel
C=func2mat(@(z) conv2(z,k,'same'), zeros([rowA,colA]));


x = optimvar("x",[rowA*colA,maxA],"Type","integer","LowerBound",0,"UpperBound",1);
y = optimvar("y",size(x),"Type","integer","LowerBound",0,"UpperBound",1);

Constraints.yUpper= y<=C*x;    %this makes y a binary map of pixels with at least 1 neighbor
Constraints.yLower= C*x/8 <= y;

Constraints.noOverlap = sum(x,2) == 1;
Constraints.minContent = sum(x)>=1;
Constraints.twoNeighbors_1 = x(:,2)<=y(:,1)+y(:,2)+3*(1-x(:,2));
Constraints.twoNeighbors_2 = C*x(:,2)+C*x(:,1) >=2*x(:,2);

Constraints.twosUpper = sum(x(:,2))<=rowA+colA; %just to keep the 2s sparse

Constraints.noBlocks =  noBlocks(x,rowA,colA);


prob = optimproblem("Constraints",Constraints,'ObjectiveSense','max');
prob.Objective=sum(x(:,2));

sol=solve(prob);

Solving problem using intlinprog.
LP:                Optimal objective value is 14.000000.                                            

Cut Generation:    Applied 116 clique cuts, 1 cover cut,                                            
                   3 zero-half cuts, and 4 Gomory cuts.                                             
                   Upper bound is 14.000000.                                                        

Heuristics:        Found 1 solution using ZI round.                                                 
                   Lower bound is 14.000000.                                                        
                   Relative gap is 0.00%.                                                          


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05.

A=round(reshape(sol.x*(1:maxA)', rowA,colA))

A = 8×6

1 2 5 5 3 1 2 2 5 1 5 3 4 4 1 2 4 5 5 5 4 4 2 2 4 4 5 5 2 1 3 5 4 5 1 2 1 4 1 2 4 2 1 3 2 2 2 2

A.*(A<=2)

ans = 8×6

1 2 0 0 0 1 2 2 0 1 0 0 0 0 1 2 0 0 0 0 0 0 2 2 0 0 0 0 2 1 0 0 0 0 1 2 1 0 1 2 0 2 1 0 2 2 2 2

function  con=noBlocks(x,rowA,colA)

    k0=zeros(3); k0([1,2,4])=1; 
    
    for i=4:-1:1
    
      k=rot90(k0,i);
      K=func2mat(@(z) conv2(z,k,'same'), zeros([rowA,colA]));
    
      con{i}= K*x(:,2)<=9*(1-x(:,2))+2*x(:,2);
    
    
    end
    con=[con{:}];

end

[b999]

Thank you very much for helping!

I will check both solutions and try to understand them. It seems like that both are working, but I am not understand every part of the solutions fully yet. Maybe I will get in touch again.

[b999]

I am suggesting, that the 9 in

 con{i}= K*x(:,2)<=9*(1-x(:,2))+2*x(:,2)
 

could be any number if it is bigger than 2?

Why do you need this part:

 - 4*(1-x(:,2)

in the constraint

 Constraints.conTwos_3 =  z >= C*x(:,2) - 4*(1-x(:,2))
 

Wouldn't that be enough:

 Constraints.conTwos_2 =  z <= C*x(:,2);
Constraints.conTwos_3 =  z >= C*x(:,2);

It does not work with this two constraints, but i don't get why I have to use exactly this expression.

[Matt J]

I am suggesting, that the 9 in .... could be any number if it is bigger than 2?

Yes.

Why do you need this part: - 4*(1-x(:,2))

Because we don't want z == C*x(:,2). We want z==x(:,2).*(C*x(:,2)). However, since that would be a nonlinear constraint, we need to look for an alternative, linear formulation.

[b999]

Hi,

I am asking myself why do you used the varibale y.

Instead of using this:

prob.Constraints.One_not_nextto_Three =  x(:,1) +  y(:,3) <=1; % 1 shouldn't be next to 3
prob.Constraints.One_not_nextto_Four =  x(:,1) +  y(:,4) <=1; % 1 shouldn't be next to 4
prob.Constraints.One_not_nextto_Six =  x(:,1) +  y(:,6) <=1; % 1 shouldn't be next to 6
prob.Constraints.One_not_nextto_Seven =  x(:,1) +  y(:,7) <=1; % 1 shouldn't be next to 7
prob.Constraints.One_not_nextto_Eight =  x(:,1) +  y(:,8) <=1; % 1 shouldn't be next to 8
prob.Constraints.One_not_nextto_One =  x(:,1) +  y(:,1) <=1; % 1 shouldn't be next to 1

why didn't you use this expressons ( yo with C*x instead of):

prob.Constraints.One_not_nextto_Three =  x(:,1) +   C*x(:,3) <=1; % 1 shouldn't be next to 3
prob.Constraints.One_not_nextto_Four =  x(:,1) +  C*x(:,4) <=1; % 1 shouldn't be next to 4prob.Constraints.One_not_nextto_Six =  x(:,1) +   C*x(:,6) <=1; % 1 shouldn't be next to 6
prob.Constraints.One_not_nextto_Seven =  x(:,1) +   C*x(:,7) <=1; % 1 shouldn't be next to 7
prob.Constraints.One_not_nextto_Eight =  x(:,1) +   C*x(:,8) <=1; % 1 shouldn't be next to 8
prob.Constraints.One_not_nextto_One =  x(:,1) +   C*x(:,1) <=1; % 1 shouldn't be next to 1

I think both works somehow, but the soltions are different. Why is it so? Is there a reason for choosing y when you want the constraint that one integer shouldn't be next to another?

[Matt J]

prob.Constraints.One_not_nextto_Three = x(:,1) + C*x(:,3) <=1

When x(:,1)=0, this will force C*x(:,3) to be <=1. In other words, all pixels A(i,j) without a 1 will now be required to have one or fewer neighboring 3s. That would be an additional restriction that you didn't have before, and may not want.

However, you could eliminate the need for y by rewriting as,

prob.Constraints.One_not_nextto_Three =  4*x(:,1) +  C*x(:,3) <=4

This will eliminate neighboring 3s when x(:,1)=1. Additionally, though, when x(:,1)=0, it only implies that C*x(:,3)<=4, which is no additional restriction at all because pixels are always fundamentally limited to at most four neighbors.

[b999]

Okay yes, it makes sense. thank you.

1) I thought about the constraints or objective function, that specific integers should have a certain distance to each other or the distance of two integers should be maximized/minimized.

For the Constraint I thought that one can use

prob.Constraints.One_not_nextto_Three =  x(:,1) +  y(:,3) <=1; 
prob.Constraints.One_not_nextto_Three =  x(:,1) +  yy(:,3) <=1;
prob.Constraints.One_not_nextto_Three =  x(:,1) +  yyy(:,3) <=1;

that the 3 should not be the neighbour (y) and not be the neighbour of the neighbour (yy) of 1 and so on (yyy, yyyy ...). As far as you want. I don't know if it is the most efficient way to implement something like that.

But I dont know how to implement an objective function which says, that the "distance" between 1 and 3 should be maximized. I think therefor you need the difference of the cells of 1 and 3...

2) And a small question about the constraints regarding the continous 2s

Constraints.twoNeighbors_1 = x(:,2)<=y(:,1)+y(:,2)+ 3*(1-x(:,2));

I understand the part, that you want the 2s to have at least one neighbour which is a 1 or a 2, so this part of the constraint:

x(:,2)<=y(:,1)+y(:,2)

but why is there this part: + 3*(1-x(:,2)). Why dou you have to set an upper limit for the parts of x(:,2) which are zero (1-x(:,2))? Couldn't you leave out this part?

[Matt J]

2) And a small question about the constraints regarding the continous 2s

This wasn't for continuous 2s. That was for the "multiple strands of 2s" problem.

[Matt J]

Couldn't you leave out this part?

I think so.

[b999]

Hey,

I thought about the two constraints for the 2s, so the "multiple strands of 2s" or the "continous 2s".

If you implement the "multiple strands of 2s" like this

Constraints.twoNeighbors_1 = x(:,2)<=y(:,1)+y(:,2)+3*(1-x(:,2));
Constraints.twoNeighbors_2 = C*x(:,2)+C*x(:,1) >=2*x(:,2);
Constraints.twosUpper = sum(x(:,2))<=rowA+colA; %just to keep the 2s sparse
Constraints.noBlocks =  noBlocks(x,rowA,colA);

and the "continous 2s" like this:

z=optimvar('z',rowA*colA,'Type','integer','Lower',0, 'Upper',2);
Constraints.conTwos_1 =  z <= 2*(x(:,2));
Constraints.conTwos_2 =  z <= C*x(:,2);
Constraints.conTwos_3 =  z >= C*x(:,2) - 4*(1-x(:,2));
Constraints.conTwos_4 =  sum(z)==2*sum(x(:,2))-2;
Constraints.conTwos_5 =  noBlocks(x,rowA,colA);

in both cases there will be a string of twos. So something like this (focussing on the 1s and 2s):

A= [0 1 0 0 0;
    0 2 2 0 0;
    0 0 2 0 0;
    0 0 2 0 0;
    0 0 2 0 0] 

In my case it would be nice, if the 2s could have branches or intersections like:

A1= [0 1 0 0 0;
     0 2 2 0 0;
     0 0 2 0 0;
     0 2 2 2 0;
     0 0 2 0 0]

or

A2= [0 1 0 0 0;
     0 2 2 0 0;
     0 0 2 0 0;
     0 0 2 2 1;
     0 0 2 0 0]

I tried to implemnt it in the constraints of "multiple strands of 2s"

Constraints.twoNeighbors_1 = x(:,2) <= y(:,1) + y(:,2) + 3*(1-x(:,2)); % 2 needs one neighbour which is 1 or 2
Constraints.twoNeighbors_2 = 4*C*x(:,2) + C*x(:,1) >= 4*x(:,2); % 2 should have max. 4 neighbours
Constraints.twosUpper = sum(x(:,2))<=rowA+colA; %just to keep the 2s sparse
Constraints.noBlocks =  noBlocks(x,rowA,colA);

But this doesn't work...i thionk there is maybe a simple solution for this in both cases.

One Intersection could also lead into a 1 (like A2), but it does not have to necessarily.

[Matt J]

In my case it would be nice, if the 2s could have branches or intersections like:

I don't think there is anything preventing branches. As you can see back here, we had a branch

[b999]

But this branch could only happen because the 2 (which is the "branch") is next to a 1. I want that there could be a branch anaywhere, not only when the 2 has a 1 as a neighbour

For example like this:

A1= [0 1 0 0 0;
     0 2 2 3 0;
     0 3 2 3 0;
     4 2 2 2 3;
     0 0 2 0 0]

I don't think there is anything preventing branches.

I don't want to prevent them, I want that they could happen along the 2s, like a road with branches or intersections.

[Matt J]

I see. That does seem difficult.

[b999]

Yes, i thought of adding something to this Constraint ( + ...)

Constraints.twoNeighbors_2 = C*x(:,2)+C*x(:,1) + ... >=2*x(:,2);
    

But I don't konw if it is the right approach. When I am trying this I often loose the connection to a 1. Maybe do you know a good approach? @Matt J?

[Matt J]

I don't, I'm afraid. It seems much harder when you allow for branches.

[b999]

Hi @Matt J,

I have a question regarding the function noBlocks. The constraint is:

con{i}= K*x(:,2)<=9*(1-x(:,2))+2*x(:,2);

The part

K*x(:,2)<= 2*x(:,2);

says, that every 2 should have maximum two 2s as neighbours around. (That's how I understand it)

But why do you need this part:

9*(1-x(:,2))

It defines the cells where is no 2. I don't understand why you implement this part. Thank you.

[Matt J]

So that no bound is placed on locations where there are no 2s.