3
$\begingroup$

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.

$\endgroup$
5
  • 2
    $\begingroup$ I think you can just use Reduce, but it's hard to say without knowing what the parameters are. $\endgroup$ Commented Aug 13, 2019 at 20:57
  • $\begingroup$ 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. $\endgroup$ Commented Aug 13, 2019 at 21:00
  • $\begingroup$ Question updated in directions to your comments. $\endgroup$ Commented Aug 13, 2019 at 22:13
  • 1
    $\begingroup$ Your e matrix has dimensions $2 \times 34$ and you're multiplying by a vector of length 32. You should fix that. $\endgroup$ Commented Aug 13, 2019 at 22:48
  • $\begingroup$ It is fixed now and the region contains at least one solution. $\endgroup$ Commented Aug 14, 2019 at 7:47

1 Answer 1

6
$\begingroup$

Read the matrices:

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[c.X[[;; 32]] <= d], 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] 

enter image description here

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}] 

enter image description here

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ Commented Aug 14, 2019 at 7:48
  • $\begingroup$ 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. $\endgroup$ Commented Aug 14, 2019 at 8:36
  • $\begingroup$ 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$. $\endgroup$ Commented Aug 14, 2019 at 9:03
  • $\begingroup$ 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. $\endgroup$ Commented Aug 14, 2019 at 9:44
  • $\begingroup$ 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. $\endgroup$ Commented Aug 14, 2019 at 10:02

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.