Why does this code return different results when use different programs?

Question

This code objective value 568 in python gekko but this code's MATLAB version gives 70. I don't understand why. Maybe about with solver options?

from gekko import GEKKO    
import numpy as np

m = GEKKO(remote=False)
capacity=np.array([70, 55, 51, 43, 41, 80])
demand=np.array([60, 57, 62, 38, 70])
cost=np.array([[5, 4, 5, 7, 2],
   [2, 9, 2, 6, 3],
   [6, 5, 1, 7, 9],
    [7, 3, 9, 3, 5],
    [4, 8, 7, 9, 7],
    [2, 5, 4, 2, 1]])

x = m.Array(m.Var,(6,5),lb=0,integer=True)

for i in range(6):
     for j in range(5):
           m.Minimize(cost[i,j]*x[i,j])

for i in range(6):
    m.Equation(m.sum(x[i,:])<=capacity[i])    


for j in range(5):
    m.Equation(m.sum(x[:,j])==demand[j])     


m.options.solver = 1
m.solve()
print('Objective Function: ' + str(m.options.objfcnval))
print(x)

Here is MATLAB version of code that give fval=70. Maybe I wrote the code wrong, but it seems mostly same, I cant understand why. Thanks for help.

clear
clc
capacity_points=6;
demand_points=5;
capacity=[70 55 51 43 41 80];
demand=[60 57 62 38 70];
cost=[5 4 5 7 2;
    2 9 2 6 3;
    6 5 1 7 9;
    7 3 9 3 5;
    4 8 7 9 7;
    2 5 4 2 1];
x=optimvar('x',capacity_points,demand_points,'Type','integer','LowerBound',0);
%expr=optimexpr;
expr=0;
for i=1:capacity_points
    for j= 1:demand_points
        expr=expr+cost(i,j)*x(i,j);
    end
end

%const1=optimconstr(capacity_points);
for i=1:capacity_points
    const1=sum(x(i,:)) <=capacity(i);
end
%const2=optimconstr(demand_points);
for j=1:demand_points
    const2=sum(x(:,j)) ==demand(j);
end

prob=optimproblem;
prob.Objective=expr;
prob.Constraints.const1=const1;
prob.Constraints.const2=const2;
sol=solve(prob)
[sol,fval] = solve(prob)

Answer 1

It appears that the solution in Gekko is correct. You may want to post the MATLAB code as well so that we can see if there are any differences between the problem statements.

Here are some general troubleshooting strategies when you want to determine if a solution is optimal:

• Determine if the problem is non-convex. This problem is a linear programming (LP) problem that fits into the form minimize c x subject to A x = b and A x < b. LPs are convex so a local solution is also the global solution. If the problem is non-convex then you can try different initial guesses to see if the solution changes or else use a global optimizer that automates that for you.

• Try different solvers. You can set up a loop and change the solver with m.options.SOLVER. If you don't want to change your code then use m.options.SOLVER = 0 to try all of the solvers. Here is the output for your problem with SOLVER=0.

 Solver           Objective  Solution Time  Status
 -------------- ------------ ------------- ---------
 APOPT (v1.0)    5.68000E+02         0.036 Success  
 BPOPT (v1.0)    5.68001E+02         0.016 Success  
 IPOPT (v3.12)   5.68000E+02         0.017 Success  
 IPOPT (v2.3)    0.00000E+00         0.000 Skip     
 SNOPT (v6.1)    0.00000E+00         0.000 Skip     
 MINOS (v5.5)    0.00000E+00         0.000 Skip     
 -------------- ------------------------------------

The freely available solvers all came to the same solution so this is an indication that multiple solver methods all achieve the same consensus.

• For anything but small-scale problems (<100 variables), it can be difficult to understand the optimization solution. Another strategy is to reduce your problem complexity until you can derive an analytic solution or until the solution is easy to verify. For your case, the columns must equal the demand while the row must stay under the capacity constraint.

capacity=np.array([70, 55, 51, 43, 41, 80])
demand  =np.array([60, 57, 62, 38, 70])
cost    =np.array([[5, 4, 5, 7, 2],
                   [2, 9, 2, 6, 3],
                   [6, 5, 1, 7, 9],
                   [7, 3, 9, 3, 5],
                   [4, 8, 7, 9, 7],
                   [2, 5, 4, 2, 1]])

If I pick out the minimum cost item for each column to meet the demand with np.min(cost,axis=0), it is [2 3 1 2 1]. If I multiply that by the demand with np.dot(np.min(cost,axis=0),demand) it gives an objective function value of 499. This is the minimum objective function without the capacity constraint. It shows that an objective of 70 is not possible for this problem unless the demand or cost is reduced.

Although these strategies are specific to your problem, they can be applied to other optimization problems to diagnose and troubleshoot a non-intuitive solution. If the solver says that it has found a solution then the Karush-Kuhn-Tucker conditions are satisfied for optimality and it is at least a local solution.

Response to Edit

The problem with the MATLAB code is that only the last constraint is enforced because const1 and const2 are redefined every loop. This gives the solution sol.x.

     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0    70

You need to include all of the equality and inequality constraints in MATLAB.

clear
clc
capacity_points=6;
demand_points=5;
capacity=[70 55 51 43 41 80];
demand=[60 57 62 38 70];
cost=[5 4 5 7 2;
    2 9 2 6 3;
    6 5 1 7 9;
    7 3 9 3 5;
    4 8 7 9 7;
    2 5 4 2 1];
x=optimvar('x',capacity_points,demand_points,'Type','integer','LowerBound',0);
%expr=optimexpr;
expr=0;
for i=1:capacity_points
    for j= 1:demand_points
        expr=expr+cost(i,j)*x(i,j);
    end
end

%const1=optimconstr(capacity_points);
for i=1:capacity_points
    const1(i)=sum(x(i,:)) <=capacity(i);
end
%const2=optimconstr(demand_points);
for j=1:demand_points
    const2(j)=sum(x(:,j)) ==demand(j);
end

prob=optimproblem;
prob.Objective=expr;
prob.Constraints.const1 = const1;
prob.Constraints.const2 = const2;
sol=solve(prob)
[sol,fval] = solve(prob)

This gives the same objective function value of 568 as gekko.