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. – Carl Woll Aug 13 '19 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. – Edmund Aug 13 '19 at 21:00
• Question updated in directions to your comments. – TEX Aug 13 '19 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. – Carl Woll Aug 13 '19 at 22:48
• It is fixed now and the region contains at least one solution. – TEX Aug 14 '19 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. – TEX Aug 14 '19 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. – Roman Aug 14 '19 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$. – TEX Aug 14 '19 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. – Roman Aug 14 '19 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. – TEX Aug 14 '19 at 10:02