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

Question

b999 2023-9-15

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2021461-problem-based-optimization-constraint-sum-x-1-1-doesn-t-work

评论： Matt J 2023-10-29

采纳的回答： Matt J

在 MATLAB Online 中打开

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?

2 个评论
显示无隐藏无

Matt J 2023-9-15

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 2023-9-15

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2023-9-15

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2021461-problem-based-optimization-constraint-sum-x-1-1-doesn-t-work#answer_1311161

编辑：Matt J 2023-9-15

在 MATLAB Online 中打开

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);

76 个评论
显示 74更早的评论隐藏 74更早的评论

Matt J 2023-9-15

编辑：Matt J 2023-9-15

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 2023-9-15

编辑：b999 2023-9-15

在 MATLAB Online 中打开

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 2023-9-15

编辑：Matt J 2023-9-15

在 MATLAB Online 中打开

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 2023-9-15

编辑：b999 2023-10-18

在 MATLAB Online 中打开

@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

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 2023-9-15

在 MATLAB Online 中打开

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 2023-9-16

编辑：b999 2023-9-17

在 MATLAB Online 中打开

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;
2 2 7;
3 2 8;
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 2023-9-17

编辑：Matt J 2023-9-17

在 MATLAB Online 中打开

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 2023-9-17

在 MATLAB Online 中打开

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.

b999 2023-9-18

在 MATLAB Online 中打开

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.

b999 2023-9-19

编辑：b999 2023-9-19

在 MATLAB Online 中打开

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.

b999 2023-9-19

编辑：b999 2023-9-19

在 MATLAB Online 中打开

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));

b999 2023-9-20

编辑：b999 2023-9-20

在 MATLAB Online 中打开

@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;
3 3 6;
4 4 2;
0 7 2;
0 5 1]
A_okay = 5×4
   2     2     8
   3     3     6
   4     4     2
   0     7     2
   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 2023-9-22

编辑：b999 2023-9-23

在 MATLAB Online 中打开

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;
 2   2   8;
 7   3   4;
 3   3   4;
 2   2   3]
A = 5×4
   2     1     5
   2     2     8
   7     3     4
   3     3     4
   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 2023-9-23

编辑：Matt J 2023-9-23

在 MATLAB Online 中打开

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 2023-9-23

在 MATLAB Online 中打开

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 2023-9-27

在 MATLAB Online 中打开

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?

b999 2023-9-27

编辑：b999 2023-9-27

在 MATLAB Online 中打开

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?

b999 2023-9-30

编辑：b999 2023-10-3

在 MATLAB Online 中打开

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;
2 2 0 0;
0 2 0 0;
0 2 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;
2 2 0 0;
0 2 0 0;
2 2 2 0;
0 2 0 0]

or

A2= [0 1 0 0 0;
2 2 0 0;
0 2 0 0;
0 2 2 1;
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.