Optimization of a matrix with integers and constraints for a car park

4 visualizaciones (últimos 30 días)
b999
b999 el 14 de Sept. de 2023
Comentada: Bruno Luong el 18 de Sept. de 2023
Hello,
I have a question about an optimization of a matrix.
I want to simulate with the matrix A a path, which a vehicle can drive. The integers have the following meaning:
1: Start/end area
2: Road
3: parking area
0: free area
A first initial matrix would be e.g.
A=
[ 0 0 1 0 0;
3 2 2 2 3;
0 0 2 0 0;
0 0 3 0 0 ]
So a start/end area, which is connected to 3 parking areas via the roads.
Now a possible objective function would be to increase the capacity of this infrastructure (matrix). The capacities result from the parking spaces and start area
Cstart = (tw/tstart) * nstart
Cpark = (tw/tpark) * npark
Where nstart and npark are the number of respective areas and represent the variables. In this case nstart = 1 and npark = 3. The total capacity would be
Call = min (Cstart, Cpark)
The objective function would be
Maximize Call.
The boundary conditions would be
- Each start area must be connected to a street
- Each parking area must be connected to a street
- There must be at least one start surface
- There must be at least one parking space
- A starting area must not be directly next to a parking area
- There must be no dead end for the streets, means that each street must be connected to two other spaces, which are not 0 (free areas).
With the objective function and the boundary conditions a matrix should then be output, which has a maximum capacity. E.g. like this:
A =
[ 3 2 1 2 3;
3 2 2 2 3;
3 2 2 2 3;
3 2 1 2 3 ]
Now, unfortunately, I don't know how to quantify these qualitative conditions. Moreover, I also don't know how, given these conditions, to derive from the objective function and the boundary conditions to the matrix A and arrange the integers.
The question now would be is there an algorithm with which one could do something like this?
  26 comentarios
Harald
Harald el 18 de Sept. de 2023
Editada: Harald el 18 de Sept. de 2023
@Bruno Luong: the problem-based approach allows for direct optimization of matrices. I find this quite helpful.
@a9v9: the least you would need is a function that returns the number of constraint violations for any input matrix. It could then be an option to use that as penalty term for the objective function rather than as an actual constraint, i.e.
max (capacity - c*[# of violations]) with a sufficiently large number c.
Depending on the size of the matrix, it could be an alternative to try all combinations in a loop. Perhaps, this is what Bruno refers to as the repeat pattern. It will of course help if you can rule out a large number of combinations to begin with.
If you combine several objectives with weights, it is still considered a single objective.
Bruno Luong
Bruno Luong el 18 de Sept. de 2023
Editada: Bruno Luong el 18 de Sept. de 2023
When I wrote "repeat pattern"I mean just repeat the 2 x 2 pattern block
[3 2;
0 1]
ans = 2×2
3 2 0 1
as I give as example for A of size( m x n) m = 4 n=4.
If m or n is odd, then truncate the last row/column then replace 1 with 3 or vice versa so that the number of 1 == number of 3.
It seems to me this simple "repeat pattern" returns the optimal solution. I don't see the need of using optimization tool for the problem stated here. But if the problem statement is different and more compliated then optimization might need to sove it.
Try with 10 x 10, I bet that the optimization would not find solution in reasonable time.
I'm curious to see if any optimization formulation can lead to finding reliable solution, and eventually beat my repeat pattern solution (or not).
m = 11;
n = 11;
A = repmat([3 2; 0 1], ceil(m/2), ceil(n/2));
A = A(1:m, 1:n);
d=floor((nnz(A==3)-nnz(A==1))/2);
i = find(A==3 & ((1:m)'==m | (1:n)==n));
A(i(1:d)) = 1
A = 11×11
3 2 3 2 3 2 3 2 3 2 3 0 1 0 1 0 1 0 1 0 1 0 3 2 3 2 3 2 3 2 3 2 3 0 1 0 1 0 1 0 1 0 1 0 3 2 3 2 3 2 3 2 3 2 3 0 1 0 1 0 1 0 1 0 1 0 3 2 3 2 3 2 3 2 3 2 3 0 1 0 1 0 1 0 1 0 1 0 3 2 3 2 3 2 3 2 3 2 3 0 1 0 1 0 1 0 1 0 1 0
nusefulplace = min(nnz(A==3),nnz(A==1))
nusefulplace = 30

Iniciar sesión para comentar.

Respuestas (1)

Bruno Luong
Bruno Luong el 18 de Sept. de 2023
Editada: Bruno Luong el 18 de Sept. de 2023
Here is formalization with intlinprog (solver base)
m = 4;
n = 5;
% Bruno repeat pattern solution
A0 = repmat([3 2; 0 1], ceil(m/2), ceil(n/2));
A0 = A0(1:m, 1:n);
d=floor((nnz(A0==3)-nnz(A0==1))/2);
i = find(A0==3 & ((1:m)'==m | (1:n)==n));
A0(i(1:d)) = 1;
b0 = A0 == reshape(0:3,1,1,4);
b0 = double(b0(:));
efficientcy0 = min([nnz(A0==3) nnz(A0==1)])
efficientcy0 = 5
% Maximize the umber of 1 and 3
z = zeros(m,n,4);
f = z; f(:,:,[2 4]) = -1;
N = numel(f);
Aeq = [];
beq = [];
lb = zeros(N,1);
ub = 1+zeros(N,1);
A1 = z; A1(:,:,2) = -1; b1 = -1;
A3 = z; A3(:,:,4) = -1; b3 = -1;
A = [A1(:)'; A3(:)'];
b = [b1; b3];
for i=1:m
for j=1:n
% There is exactly one ID by position
k = sub2ind([m n 4], i+zeros(1,4), j+zeros(1,4), 1:4);
Ak = accumarray(k(:), 1, [N,1]);
Aeq(end+1,:) = Ak;
beq(end+1) = 1;
% 4-connected neighbor
in = [i-1, i+1, i, i];
jn = [j, j, j-1, j+1];
keep = in>=1 & in<=m & jn>=1 & jn<=n;
in = in(keep);
jn = jn(keep);
% 1 must be connected to 2
k = sub2ind([m n 4], i, j, 2);
kn = sub2ind([m n 4], in, jn, 3+zeros(size(in)));
Ak = -accumarray(kn(:), 1, [N,1]) + accumarray(k(:), 1, [N,1]);
A(end+1,:) = Ak;
b(end+1) = 0;
% 3 must be connected to 2
k = sub2ind([m n 4], i, j, 4);
kn = sub2ind([m n 4], in, jn, 3+zeros(size(in)));
Ak = -accumarray(kn(:), 1, [N,1]) + accumarray(k(:), 1, [N,1]);
A(end+1,:) = Ak;
b(end+1) = 0;
% 2 must be connected to 1, 2, or 3
k = sub2ind([m n 4], i, j, 3);
kn1 = sub2ind([m n 4], in, jn, 2+zeros(size(in)));
kn2 = sub2ind([m n 4], in, jn, 3+zeros(size(in)));
kn3 = sub2ind([m n 4], in, jn, 4+zeros(size(in)));
kn = [kn1, kn2, kn3];
Ak = -accumarray(kn(:), 1, [N,1]) + accumarray(k(:), 1, [N,1]);
A(end+1,:) = Ak;
b(end+1) = 0;
% 1 must not be connected to 3
k = sub2ind([m n 4], i, j, 2);
for l=1:length(in)
kn = sub2ind([m n 4], in(l), jn(l), 4);
Ak = accumarray(kn(:), 1, [N,1]) + accumarray(k(:), 1, [N,1]);
A(end+1,:) = Ak;
b(end+1) = 1;
end
end
end
% Balance between 1 and 3
Ak = repmat(reshape([0 1 0 -1],1,1,4),m,n,1);
A(end+1,:) = +Ak(:); b(end+1) = 1;
A(end+1,:) = -Ak(:); b(end+1) = 1;
bin = intlinprog(f, 1:N, A, b(:), Aeq, beq(:), lb, ub, b0);
LP: Optimal objective value is -17.055556. Cut Generation: Applied 6 Gomory cuts, and 8 zero-half cuts. Lower bound is -15.000000. Relative gap is 45.45%. Heuristics: Found 2 solutions using diving. Upper bound is -13.000000. Relative gap is 14.29%. Cut Generation: Applied 2 strong CG cuts, and 12 zero-half cuts. Lower bound is -14.000000. Relative gap is 7.14%. Branch and Bound: nodes total num int integer relative explored time (s) solution fval gap (%) 29 0.18 3 -1.300000e+01 0.000000e+00 Optimal solution found. Intlinprog stopped 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.
bin = reshape(round(bin), [m,n,4]);
A = sum(bin.*reshape(0:3,1,1,4),3)
A = 4×5
2 3 3 2 3 3 3 2 1 2 3 2 1 1 1 3 2 1 2 1
efficientcy = min([nnz(A==3) nnz(A==1)])
efficientcy = 6
  2 comentarios
b999
b999 el 18 de Sept. de 2023
@Bruno Luong Thanks for helping me. I will check the code and try to understand it. Due to my knowledge in optimization and matlab it could take some time...i am not familiar with many operations you use in this code.
Is it possible in general to add some more integers to the Matrix, e.g. that you have 0,1,2,3,4,5,6? Or do you have to rewrite the whole code in this case?
One small question at first: why aren't you defining A2 in this part?
A1 = z; A1(:,:,2) = -1; b1 = -1;
A3 = z; A3(:,:,4) = -1; b3 = -1;
A = [A1(:)'; A3(:)'];
Bruno Luong
Bruno Luong el 18 de Sept. de 2023
Because your constraints include these
  • There must be at least one start surface
  • There must be at least one parking space
But you do not tell "There must be at least one road"

Iniciar sesión para comentar.

Categorías

Más información sobre Multiobjective Optimization en Help Center y File Exchange.

Productos


Versión

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by