Plot feasible region of a high-dimensional linear programming along some dimensions
4 ビュー (過去 30 日間)
古いコメントを表示
Let
be a
vector of unknowns. Consider the linear system
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233837/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233838/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233839/image.png)
where
are matrices and vectors of known parameters with appropriate dimensions. They are attached in .mat format in the last comment to the question. Let
.
is bounded.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233840/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233841/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233793/image.png)
My objective: I would like to plot the region of values of
for which there exists
such that
satisfies the linear system above. Let us call such a region
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233842/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233843/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233844/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233797/image.png)
I initially thought about the following strategy to reach my objective.
1) Find the extreme points of
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233793/image.png)
2) Transform the extreme points of
using E to get a description of the extreme points of
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233793/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233797/image.png)
3) Use fill to colour the area traced by the extreme point of
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233797/image.png)
[V,nr,nre]=lcon2vert(C,d,A,b,[])
which stops with a long error message that I'm struggling to intepret. On top of that, even if I was able to make the function lcon2vert running, I doubt that 32 unknowns is a manageable size.
I would like to know your opinion on this. Could you think of other strategies? For example an algorithm that find random points in
across each iteration may be an idea: as the number of iteration goes to infinity, the scatter of the points gets closer and closer to
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233797/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233797/image.png)
Clarification: The linear system above has at least one solution (see possible_solution_complete.mat)
7 件のコメント
採用された回答
Bruno Luong
2019 年 8 月 14 日
編集済み: Bruno Luong
2019 年 8 月 14 日
Something is wrong in my implementation of understanding of your problem, but there is no feasible point found by my code.
A=load('A.txt');
b=load('b.txt');
C=load('C.txt');
d=load('d.txt');
E=load('E.txt');
Aeq = [A, zeros(size(A,1),2); ...
E, -eye(2)
];
beq = [b; ...
zeros(size(E,1),1)
];
N = 360;
Rb = nan(2,N);
for i=0:N-1
theta = 2*pi/N*i;
fi = [zeros(32,1); cos(theta); sin(theta)];
[xi,~, exitflag] = linprog(fi, ...
[C, zeros(size(C,1),2)], d, ...
Aeq, beq, [], []);
if exitflag>=0
Rb(:,i+1) = xi([33 34]);
end
end
Edit: this code no longer applicable, see updated code using mat files below
18 件のコメント
Matt J
2019 年 8 月 16 日
It indicates the boundary/domain is a line. But again there might be just an artifact of discretization (or not). Who knows. I increases N to 3600, I still get a line.
This can also be seen by eliminating the equality constraints A*y=b. The solutions have the form y=Na*c+y0 where Na=null(A) and y0 is the given feasible point. The solutions z=[x33,x34] then have the form z=E*(Na*c+y0). However,
>> Na=null(A); rank(E*Na,1e-12)
ans =
1
So, z lies in some 1-dimensional space.
その他の回答 (2 件)
Bruno Luong
2019 年 8 月 14 日
編集済み: Bruno Luong
2019 年 8 月 14 日
What is your motivation to do that?
You clearly underestimate such task, simply because the number of vertices of simplex grow very fast with the dimension.
I see your inequality constrainst are 72 (C matrix) assuming they are independent and number of variable is 34. Then the number of vertices can go as big as
>> nchoosek(72,34)
Warning: Result may not be exact. Coefficient is greater than 9.007199e+15 and is only accurate to 15 digits
> In nchoosek (line 92)
ans =
3.9656e+20
If you project on subspace defined by A ane E (20+2 rows) then you get a subspace of dimension 12, then the number of vertices is still very large
>> nchoosek(72,34-22)
ans =
1.5363e+13
To store all of them you need
nchoosek(72,34-22)*8/1e12
~123 Tb.
Each of the vertices require to solve a linear system of 34 x 34.
There is very little chance that you could find a method to do it in reasonably time and with enough accuracy.
Generate a random points to cover such shape? Good luck!!!
2 件のコメント
Bruno Luong
2019 年 8 月 14 日
編集済み: Bruno Luong
2019 年 8 月 14 日
R is projection of P in 2-dimensional subspace, so I think there is no short-cut,
Now if you want an approximation of R you can discretize:
f_i := [0,0,...,0, cos(theta_i), sin(theta_i)].'
where
theta_i is 2*pi*i/N; i=0,1,2....,N-1 and N large enough.
Finding the N solutions of N LP
x_i = argmin f_i'*x such that
equality/inequalities constraints with matrices A, C, E are meet.
After taking coordinated 33/34 of {x_i} you will getsome discretization of the boundary of R, but there is no warranty what so ever that descretization is the true shape.
Matt J
2019 年 8 月 17 日
I tried to use this without success ... which stops with a long error message that I'm struggling to intepret.
The error message indicates that the region described by C*y<=d is not solid, i.e., it lies in some strict subspace of R^32, to within numerical precision. This is a numerically unstable way of specifying your constraints. It means that some of the inequalities in the system C*x<=d should really be expressed as equalities.
If you think this assessment is inaccurate, try to present an interior point y, one that strictly satisifes the inequalities C*y<d.
0 件のコメント
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!