Solving system of equations with constraints

I am solving a system of linear equations with constraints using linprog.

Aeq = [-1 0 1 0 0 0 0;1 -1 0 -1 0 0 0;0 1 0 0 1 0 0;-1 0 0 0 0 1 0;1 0 0 0 0 0 1];
rank(Aeq)
size(A)
beq = [0 0 0 0 0];
size(b)
A=[];
b=[];
f = [0 0 0 -1 0 0 0];
lb = [0 -10 0 0 -10 -10 -10];
ub = [10 10 10 10 10 10 10];
x = linprog(f,A,b,Aeq,beq,lb,ub)

The optimal solution that I obtain is,

x =

     0
   -10
     0
    10
    10
     0
     0

From what I observe, the values in the solution vector are mostly the values of the bounds.

Could someone suggest how to obtain other solutions?

For instance,

x = [ 4.857 ; 5.143; 4.857;10;5.143;4.857;4.857] is also a solution.

I am not sure how other solutions can be obtained.

Answer 1

[John D'Errico]

Linear programming returns the optimal solution. You want another optimal solution? How many optimal solutions are there? :) Ok, lets get serious.

You have a very simple LP problem. Maximiize variable 4, subject to a set of equality constraints, and bound constraints.

Linprog finds this solution:

xlinprog =
     0
   -10
     0
    10
    10
     0
     0

with x(4) at the upper bound. So, what flexibility do you have?

Aeq
Aeq =
    -1     0     1     0     0     0     0
     1    -1     0    -1     0     0     0
     0     1     0     0     1     0     0
    -1     0     0     0     0     1     0
     1     0     0     0     0     0     1

Aeq*x == 0 is a system of 5 equations in 7 unknowns. So there will be two degrees of freedom in there.

format long g
Anull = null(Aeq)
Anull =
         0.297698713811669           -0.354487221322
        -0.342890338788254         -0.48945793762951
         0.297698713811669           -0.354487221322
         0.640589052599922         0.134970716307511
         0.342890338788254          0.48945793762951
         0.297698713811669           -0.354487221322
        -0.297698713811669            0.354487221322

That is, we can add any linear combination of the two columns of null(Aeq) to the solution, and get another solution that lies on the equality constraint. Efectively, we can create a family of quasi-solutions that look like

xlinprog + s*Anull(:,1) + t*Anull(:,2)

The problem is to do so, while x(4) remains unchanged. What does that mean? We can now formulate one other simple equation:

Anull(4,1)*s + Anull(4,2)*t == 0

Thus, given two parameters s and t, we can add any amount of these two vectors to the linprog solution, as long as that additional equality holds true. In turn, that lets us eliminate t, as

t = -s*Anull(4,1)/Anull(4,2)

So the family of all "solutions" to your problem are of the form

xlinprog + s*Anull(:,1) - s*Anull(:,2)*Anull(4,1)/Anull(4,2)

That should make it clear there is a one dimensional family of probably infinitely many solutions, but those other solutions will not lie at a vertex of the solution space, but be interior to that domain. Since linprog will return a solution that always lies at a boundary vertex, it gave you what it did.

xfamily = @(s) xlinprog + s*Anull(:,1) - s*Anull(:,2)*Anull(4,1)/Anull(4,2)
lb = [0 -10 0 0 -10 -10 -10]';
ub = [10 10 10 10 10 10 10]';

Now, we can quickly observe that xfamily(0) is the linprog solution, and that the entire set of validly equivalent solutions lies in the interval on s of [0,5.05014177670167].

[lb,xfamily(0), xfamily(5.05014177670167),ub]
ans =
            0            0           10           10
          -10          -10   7.1054e-15           10
            0            0           10           10
            0           10           10           10
          -10           10  -7.1054e-15           10
          -10            0           10           10
          -10            0          -10           10

Anything between those two values will result in a validly equivalent solution, that remains within the bounds.

Easy peasy.

Answer 2

[Torsten]

Linear optimization codes that use the simplex method always return corner points of the LP problem as solution.

Maybe

options = optimoptions('linprog','Algorithm','interior-point')

can help.

[Deepa Maheshvare]

Thank you, I tried

options = optimoptions('linprog','Algorithm','interior-point ') 
x = linprog(f,A,b,Aeq,beq,lb,ub,options)

But, it returns the corner point again.