You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Problem-based Optimization: Constraint sum(x==1) >=1 doesn't work.
11 views (last 30 days)
Show older comments
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 Comments
Matt J
on 15 Sep 2023
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?
Accepted Answer
Matt J
on 15 Sep 2023
Edited: Matt J
on 15 Sep 2023
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 Comments
b999
on 15 Sep 2023
Edited: b999
on 17 Sep 2023
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
on 15 Sep 2023
Edited: Matt J
on 15 Sep 2023
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
on 15 Sep 2023
Edited: b999
on 15 Sep 2023
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
on 15 Sep 2023
Edited: Matt J
on 15 Sep 2023
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
on 15 Sep 2023
Edited: b999
on 15 Sep 2023
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
on 15 Sep 2023
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
on 15 Sep 2023
Edited: Matt J
on 15 Sep 2023
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
on 15 Sep 2023
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
on 15 Sep 2023
Edited: Torsten
on 15 Sep 2023
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
on 15 Sep 2023
Edited: Matt J
on 15 Sep 2023
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
on 15 Sep 2023
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
on 15 Sep 2023
Edited: b999
on 18 Oct 2023
@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));
b999
on 15 Sep 2023
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
on 15 Sep 2023
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
on 16 Sep 2023
Edited: b999
on 17 Sep 2023
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
on 17 Sep 2023
Edited: Matt J
on 17 Sep 2023
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,
b999
on 17 Sep 2023
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
on 17 Sep 2023
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
on 18 Sep 2023
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
on 18 Sep 2023
Edited: Matt J
on 18 Sep 2023
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
on 18 Sep 2023
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
on 18 Sep 2023
Well, a guaranteed way to avoid it is to provide a feasible initial guess, assuming you know one.
b999
on 18 Sep 2023
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
on 18 Sep 2023
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
on 18 Sep 2023
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.
b999
on 19 Sep 2023
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?
b999
on 19 Sep 2023
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
on 19 Sep 2023
p.x = x0;
p.y = y0;
% + all other optimvars you have defined
inf_test = infeasibility(prob.Constraints.Eight_nextto_Two,p)
Matt J
on 19 Sep 2023
Edited: Matt J
on 19 Sep 2023
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
on 19 Sep 2023
Edited: b999
on 19 Sep 2023
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.
b999
on 19 Sep 2023
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
on 19 Sep 2023
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
on 19 Sep 2023
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
on 19 Sep 2023
Edited: b999
on 19 Sep 2023
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
on 19 Sep 2023
Edited: b999
on 19 Sep 2023
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
on 19 Sep 2023
Well, understanding why something is infeasible requires that we see both A and the code you used to generate sol0.
b999
on 19 Sep 2023
Edited: b999
on 19 Sep 2023
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
on 19 Sep 2023
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
on 20 Sep 2023
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.
Matt J
on 20 Sep 2023
Edited: Matt J
on 20 Sep 2023
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
on 20 Sep 2023
Edited: b999
on 20 Sep 2023
@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
on 22 Sep 2023
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
on 22 Sep 2023
Edited: Matt J
on 22 Sep 2023
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
on 22 Sep 2023
Edited: Matt J
on 22 Sep 2023
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
on 22 Sep 2023
Edited: b999
on 23 Sep 2023
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
on 23 Sep 2023
Edited: Matt J
on 23 Sep 2023
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
on 23 Sep 2023
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
on 24 Sep 2023
Edited: b999
on 24 Sep 2023
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
on 24 Sep 2023
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
on 27 Sep 2023
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
on 27 Sep 2023
Edited: Matt J
on 27 Sep 2023
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
on 27 Sep 2023
Edited: b999
on 27 Sep 2023
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
on 30 Sep 2023
Edited: b999
on 3 Oct 2023
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
on 3 Oct 2023
In my case it would be nice, if the 2s could have branches or intersections like:
b999
on 3 Oct 2023
Edited: b999
on 3 Oct 2023
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.
b999
on 3 Oct 2023
Edited: b999
on 5 Oct 2023
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?
b999
on 29 Oct 2023
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.
More Answers (0)
See Also
Categories
Find more on Get Started with Optimization Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)