# Plot feasible region of a high-dimensional LP only along some dimensions

I have a question on the possibility of using Mathematica to plot a convex closed region satisfying a linear system of equalities and inequalities.

Let me first present the linear system. Let $$x\equiv (x_1,...,x_{34})$$ be a $$34\times 1$$ vector of unknowns. Let $$A,C,E,b,d$$ be matrices of known parameters with appropriate dimensions. The linear system is

$$1) A\times (x_1,...,x_{32})=b$$ ,

$$2)C\times(x_1,...,x_{32}) \leq d$$,

$$3)E\times(x_1,...,x_{32})-(x_{33},x_{34})=0$$,

The matrices are here https://filebin.net/e7k3749uxd2f1dg4, in .mat format. They are too big to be reported.

My objective: I would like to plot the region of values of $$(x_{33},x_{34})$$ for which there exists $$(x_1,...,x_{32})$$ such that $$(x_1,...,x_{34})$$ satisfies $$1),2),3)$$.

Question: Can Mathematica allow me to plot the desired $$2$$-D region? The tricky part here is to explore the set of solution of $$1),2)$$ with respect to $$(x_1,...,x_{32})$$. I typically use Matlab which, however, to the best of my knowledge, does not have packages doing what I want due to the high dimension of the problem.

Clarification: I have never used Mathematica (hence, I don't have a code of attempts to show you), but I'd be happy to start studying it if you tell me that it can help me with my question.

Clarification 2: there exists at least one value of $$(x_1,...,x_{34})$$ satisfying 1),2),3). It is here https://filebin.net/e7k3749uxd2f1dg4 under the name possible_solution_complete.mat.

• I think you can just use Reduce, but it's hard to say without knowing what the parameters are. Commented Aug 13, 2019 at 20:57
• A good start would be to provide the set matrices and constraints. Also, please include code of what you have tried and what is not working for you. Commented Aug 13, 2019 at 21:00
– Star
Commented Aug 13, 2019 at 22:13
• Your e matrix has dimensions $2 \times 34$ and you're multiplying by a vector of length 32. You should fix that. Commented Aug 13, 2019 at 22:48
• It is fixed now and the region contains at least one solution.
– Star
Commented Aug 14, 2019 at 7:47

a = Import["A.mat"][[1]] // Round;
c = Import["C.mat"][[1]] // Round;
e = Import["E.mat"][[1]];
b = Import["b.mat"] // Flatten;
d = Import["d.mat"] // Flatten // Round;


construct the system of equalities and inequalities:

X = Array[x, 34];
S = Join[Thread[a.X[[;; 32]] == b],
Thread[e.X[[;; 32]] - {X[[33]], X[[34]]} == 0]];


Experimental mathematics: find 100 random instances that satisfy S:

j = FindInstance[S, X, Reals, 10^2];


For each instance, plot $$(x_{33},x_{34})$$:

J = X[[-2 ;;]] /. j;
ListPlot[J, AspectRatio -> Automatic]


Notice that these points lie on the line $$x_{33}+x_{34}=1$$:

MinMax[Total /@ J]
(*    {0.9999999999999998, 1.0000000000000002}    *)


Find the vertices in the $$x_{33}-x_{34}$$ plane by minimizing/maximizing these parameters under the constraints:

P1 = X /. Minimize[{x[33], S}, X][[2]]
(*    {0., 0.304716, 0., 0.044609, 0., 0.544515, 0., 0.10616, 0.304716,
0., 0.044609, 0., 0., 0.544515, 0.10616, 0., 0., 0.304716, 0.,
0.044609, 0., 0.544515, 0., 0.10616, 0., 0.304716, 0., 0.044609,
0., 0.544515, 0., 0.10616, 0.0505921, 0.949408}    *)

P2 = X /. Maximize[{x[33], S}, X][[2]]
(*    {0.304716, 0., 0.044609, 0., 0.044609, 0.499906, 0.10616, 0.,
0.304716, 0., 0.044609, 0., 0.544515, 0., 0.10616, 0., 0.,
0.304716, 0.044609, 0., 0., 0.544515, 0., 0.10616, 0.304716, 0.,
0.044609, 0., 0.044609, 0.499906, 0.10616, 0., 0.415966, 0.584034}    *)

P3 = X /. Minimize[{x[34], S}, X][[2]]
(*    {0.304716, 0., 0.044609, 0., 0.044609, 0.499906, 0.10616, 0.,
0.304716, 0., 0.044609, 0., 0.544515, 0., 0.10616, 0., 0.,
0.304716, 0.044609, 0., 0., 0.544515, 0., 0.10616, 0.304716, 0.,
0.044609, 0., 0.044609, 0.499906, 0.10616, 0., 0.415966, 0.584034}    *)

P4 = X /. Maximize[{x[34], S}, X][[2]]
(*    {0., 0.304716, 0., 0.044609, 0., 0.544515, 0., 0.10616, 0.304716,
0., 0.044609, 0., 0., 0.544515, 0.10616, 0., 0., 0.304716, 0.,
0.044609, 0., 0.544515, 0., 0.10616, 0., 0.304716, 0., 0.044609,
0., 0.544515, 0., 0.10616, 0.0505921, 0.949408}    *)


There are actually only two inequivalent points:

P1 == P4
(*    True    *)

P2 == P3
(*    True    *)


Therefore, the domain you're looking for is the line from P1 to P2, projected onto the $$x_{33}-x_{34}$$ plane:

{P1[[-2 ;;]], P2[[-2 ;;]]}
(*    {{0.0505921, 0.949408}, {0.415966, 0.584034}}    *)


You can plot the domain (in general) with

ConvexHullMesh[{P1, P2, P3, P4}[[All, -2 ;;]], Frame -> True, FrameLabel -> {x33, x34}]


• Thanks, I'm sorry but I did a mistake when building the matrices! I had to convert them from Matlab format and this confused me! They are fixed now and the program has at least one solution.
– Star
Commented Aug 14, 2019 at 7:48
• Can you tell us what that one known solution is? Because I just tried with your new coefficients and still don't get any solutions. Commented Aug 14, 2019 at 8:36
• I have uploaded it under the same link. I know little of Mathematica coding but I think that maybe Thread[e.X == 0] does not correspond to $E\times (x_1,...,x_{32})-(x_{33},x_{34})=0$.
– Star
Commented Aug 14, 2019 at 9:03
• Yes I had modified it to Thread[e.X[[;; 32]] - {X[[33]], X[[34]]} == 0] before trying; still no luck. But you're right that your solution seems to work. Commented Aug 14, 2019 at 9:44
• I understood the problem I think: the .txt extension was producing wrong numbers (I haven't got yet why). Now, the files are in the .mat extension and should work fine.
– Star
Commented Aug 14, 2019 at 10:02