2
$\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$
  • 2
    $\begingroup$ I think you can just use Reduce, but it's hard to say without knowing what the parameters are. $\endgroup$ – Carl Woll Aug 13 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$ – Edmund Aug 13 at 21:00
  • $\begingroup$ Question updated in directions to your comments. $\endgroup$ – user3285148 Aug 13 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$ – Carl Woll Aug 13 at 22:48
  • $\begingroup$ It is fixed now and the region contains at least one solution. $\endgroup$ – user3285148 Aug 14 at 7:47
5
$\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$
  • $\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$ – user3285148 Aug 14 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$ – Roman Aug 14 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$ – user3285148 Aug 14 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$ – Roman Aug 14 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$ – user3285148 Aug 14 at 10:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.