# Plot feasible region of a high-dimensional linear programming along some dimensions

10 views (last 30 days)
CT on 13 Aug 2019
Answered: Matt J on 17 Aug 2019
Let be a vector of unknowns. Consider the linear system
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.
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 .
I initially thought about the following strategy to reach my objective.
1) Find the extreme points of .
2) Transform the extreme points of using E to get a description of the extreme points of .
3) Use fill to colour the area traced by the extreme point of .
I'm stuck at 1). I tried to use this without success, by running
[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 .
Clarification: The linear system above has at least one solution (see possible_solution_complete.mat)
Matt J on 14 Aug 2019
For the convenience of the community, I have reposted CT's data in a single .mat file

Bruno Luong on 14 Aug 2019
Edited: Bruno Luong on 14 Aug 2019
Something is wrong in my implementation of understanding of your problem, but there is no feasible point found by my code.
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
Matt J on 16 Aug 2019
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.

Bruno Luong on 14 Aug 2019
Edited: Bruno Luong on 14 Aug 2019
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 CommentsShowHide 1 older comment
Bruno Luong on 14 Aug 2019
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 on 17 Aug 2019
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.