How to solve the following equations using simplex method?

Question

Software Engineer here,

I am trying to find an algorithm to solve the following problem, basically I have 3 equations that you can see bellow, and all values of X, Y, Z, and Xi, Yi, Zi's are known. The only unknowns are C values that I am trying to find.

I understand Simplex Method has to be used there (or if anything else please suggest).

But I am new to simplex method, and really confused about many things, like for example what is my objective function? I understand all equalities should be changed to 2 inequalities, so that way I have 6 equations, and this can be considered my restrictions? in that case still confused about my objective function. What should I maximize or minimize if I am just trying to find a value?

If anyone can help me understand this better so I can eventually understand how to make a Tableu and solve this using a programming language, will be great. (Links to a good reads are appreciated as well, so far tried wikipedia, wasn't a good help)

am I even on the right path?

Anyway, here are the equations:

Edit: Forgot to add, all variables are between 0 and 1, which is a major constraint I guess.

Edit 2: I made some progress since yesterday, and tried to implement the simplex the way I see. (See the Tableau)(I tried to maximize for the SUM of C's as a goal)

And it kind of works! As in, it did calculate most of the cases exactly right.

Here is how I test if it was correct - I take my numbers feed to simplex, get C's, then I multiply vector of C's with the matrix back again, expecting to get the same X Y Z values I started with. If it's the same, then it worked.

Problem is, there are weird edge cases! That I can't seem to be able to wrap my brain around.

For example this values work perfectly: X = 0.06837372481822968 Y = 0.13674744963645935 Z = 0.022791240364313126

But, this values (literally almost the same) X = 0.06716471165418625 Y = 0.1343294233083725 Z = 0.022388236597180367

fail!, And fail means the resulting C values from simplex are HUGELY different (missing mid part, middle of C's are zeroes), and this when multiplied back with matrix produces different results from initial values.

How can that be? does it mean that simplex fails due to some wrong constraints or? How do I look at this?

To explain this better, take a look how resulting answer of simplex, just collapses with this little number chance (I checked and during process at some point just different pivot is chosen) Check how third line just dropped in the middle, compared to other two.

This pic kind of suggests that issue is because solutions dip under 0, for whatever reason? not sure why and how to prevent that.

Answer 1

Linear programming is about optimizing a linear objective function over a polyhedral. Ask yourself, what are the desirable property that you intend to optimize. If your goal is just to obtain a feasible point, you can just optimize $0$, the trivial objective function. Or you could have optimize the sum of the variable (which is what you have chosen).

Now, about your weird cases. The most alarming information to me is

and this when multiplied back with matrix produces different results from initial values.

If I understand this correctly and you mean $Ax =b$ no longer hold approximately, then there is some bug.

Simplex algorithm should preserve the equality constraint at all time theoretically. In practice, there could be some minor numerical differences of course. If this is not the case, then something has gone wrong. You might want to capture the first moment this condition is violated and also check your checking procedure.

I have called the simplex function from Python Scipy library to check if the reported behavior can be reproduced. Do check if I have constructed your constraint matrix correctly.

import numpy as np
from scipy import optimize
sub_A = np.array([[0.003, 0.168, 0.098, 0.122, 0.502, 0.454, 0.072, 0.003, 0, 0],[0, 0.028, 0.169, 0.503, 0.539, 0.231, 0.029, 0.001, 0 , 0], [0.015, 0.854, 0.698, 0.042, 0, 0,0,0,0,0]])
top = np.concatenate((sub_A, np.zeros((3,10))), axis = 1 )
bottom = np.concatenate((np.identity(10), np.identity(10)), axis =1)
fixed_A = np.concatenate((top, bottom), axis = 0)
fixed_C = np.concatenate((-np.ones(10), np.zeros(10)), axis = 0)

def linprogcheck(x, y, z):
    b = np.concatenate((np.array([x,y,z]), np.ones(10)), axis = 0)
    ans = optimize.linprog(fixed_C, method = 'simplex', A_eq = fixed_A, b_eq = b)
    print("check error size:")
    print(np.linalg.norm(np.matmul(fixed_A, ans.x)-b))
    return ans

ans1 = linprogcheck(0.06837372481822968 ,0.13674744963645935 , 0.022791240364313126 )
print(ans1)

ans2 = linprogcheck(0.06716471165418625, 0.1343294233083725, 0.022388236597180367)
print(ans2)

The output that I obtained shows that the equality constraint holds and the two solutions are closed to each other.

check error size:
1.4304896245381992e-17
     fun: -4.542226973252382
 message: 'Optimization terminated successfully.'
     nit: 20
   slack: array([], dtype=float64)
  status: 0
 success: True
       x: array([0.83814691, 0.        , 0.        , 0.2433104 , 0.        ,
       0.        , 0.46076966, 1.        , 1.        , 1.        ,
       0.16185309, 1.        , 1.        , 0.7566896 , 1.        ,
       1.        , 0.53923034, 0.        , 0.        , 0.        ])
check error size:
3.122502256758253e-17
     fun: -4.514192719488299
 message: 'Optimization terminated successfully.'
     nit: 20
   slack: array([], dtype=float64)
  status: 0
 success: True
       x: array([0.82330393, 0.        , 0.        , 0.23901614, 0.        ,
       0.        , 0.45187266, 1.        , 1.        , 1.        ,
       0.17669607, 1.        , 1.        , 0.76098386, 1.        ,
       1.        , 0.54812734, 0.        , 0.        , 0.        ])

Answer 2

The main reason for the jumps between unknowns $C_i$ in the simplex method is that both the goal function and the constraints of task are symmetric over the unknowns $C_i$.

Taking into account that the basic constraints can be presented as the equations of hyperplanes via axes intersections, the upper bounds of the coordinates can be calculated by the formula $$C_i^{max}= \min\left(1, \dfrac X{X_i}, \dfrac Y{Y_i}, \dfrac Z{Z_i}\right).$$ Also, can be used the goal function in the form of $$G = \sum_{i=1}^{10} (1+\varepsilon i) C_i.$$

This should improve the simplex method, both in terms of speed and in terms of sustainability.