# input variables
Constraint_Coefficients = #constraintCoefficients#;
Constraint_Constants = #constraintConstants#;
Objective_Coefficents = #objectiveCoefficients#;
Constraint_Type = #constraintTypes#;
Variables_Names = #variableNames#; # can use "x" instead of a array
Variable_Type = #variableTypes#;
Problem_Type = '#ddlMinMax#';
#define function to flip inequality sign
def flip_sign(sn):
if sn == ">=":
return "<=";
if sn == "<=":
return ">=";
#find tableau
def tableau(D,phase_variables_list,pb_type,ob_name):
# get current variables from current LP dictionary
current_variables_list=list(D.nonbasic_variables())+list(D.basic_variables());
# get current variable order
current_variables_order = list(range(len(phase_variables_list)));
for variable in phase_variables_list:
current_variables_order[current_variables_list.index(variable)] = phase_variables_list.index(variable);
# define the double array
A = [[0]*len(current_variables_list) for i in list(range(len(D.basic_variables())))];
# retrieve coefficients of current states
for i in range(len(D.basic_variables())):
for j in range(len(D.nonbasic_variables())):
A[i][current_variables_order[j]] = D.column_coefficients(current_variables_list[j])[i];
# fill the one for basic variables
for i in range(len(D.basic_variables())):
A[i][current_variables_order[len(D.nonbasic_variables())+i]] = 1;
# add solution of each row
for i in range(len(D.basic_variables())):
A[i] += [D.constant_terms()[i]];
# add 'basic variables' column
for i in range(len(D.basic_variables())):
A[i] = [D.basic_variables()[i]] + A[i];
# add objective row
objective_row =[0]*len(current_variables_list);
for j in range(len(D.nonbasic_variables())):
if pb_type == 'max':
objective_row[current_variables_order[j]] = -1*D.objective_coefficients()[j];
else:
objective_row[current_variables_order[j]] = D.objective_coefficients()[j];
if (pb_type == 'min') or (pb_type == '-max'):
objective_row += [-1*D.objective_value()];
else:
objective_row += [D.objective_value()];
obj_name = ob_name;
objective_row = [obj_name] + objective_row;
A = [objective_row] + A;
# add head row
head_row = ['basic'] + phase_variables_list + ['solution'];
A = [head_row] + A;
return A
# end of tableau function
# display tableau function
def dis_tableau(tableau):
T = tableau;
begin_code = '$$$$ \\\\begin{array}';
end_code = '\\\\hline \\\\end{array} $$$$';
col_code = '{ | c |';
for vl in T[0]:
col_code += ' l ';
col_code += '|}';
#add head row
content_code = '\\\\hline ';
for name in T[0]:
content_code += str(name) + ' & ';
content_code = content_code[:-2];
content_code += '\\\\\\\\ \\\\hline ';
for i in range(1,len(T)):
for j in range(len(T[0])):
content_code += ' ' + str(T[i][j]) ;
content_code += ' & '
content_code += ' \\\\\\\\ ';
if i == 1:
content_code += ' \\\\hline ';
tex_code = begin_code + col_code + content_code + end_code
show(html(tex_code))
#define function to make the ratios table
def ratio_table(D):
entering = D.entering();
entered_ratio = D.ratios();
A = [];
for row in entered_ratio:
row_index = list(D.basic_variables()).index(row[1]);
A += [[row[1],D.column_coefficients(entering)[row_index],D.constant_terms()[row_index],row[0]]]
A = [["Possible Leaving",str(entering) + " (entering)","Solutions","Ratio"]] + A
return A
#ratios table function end
#main
#construct LP problem class
Con_Coe = list(Constraint_Coefficients);
Con_Const = list(Constraint_Constants);
Obj_Coe = list(Objective_Coefficents);
P = InteractiveLPProblem(Con_Coe, Con_Const, Obj_Coe, Variables_Names, problem_type=Problem_Type,constraint_type=Constraint_Type,variable_type=Variable_Type,);
#check feasibility, if feasible, tell tehm to use ch1 calculator
if P.is_feasible([0]*len(Variables_Names)):
show(html('
The origin is inside the feasible set by constraints,
you may use the simplex calculator in chapter one.
'))
assert()
#retrieve base ring of the problem
Base_Ring = P.base_ring();
#retrieve variable_list
variables_list = list(P.x());
show(P)
#flip problem type from min to -max
if Problem_Type == "min":
Problem_Type = "-max"
Obj_Coe = list(-1*vector(Obj_Coe));
P = InteractiveLPProblem(Con_Coe, Con_Const, Obj_Coe, Variables_Names, problem_type=Problem_Type,constraint_type=Constraint_Type,variable_type=Variable_Type,);
# change the constraints if they are not in standard form
is_flip = false;
for i in range(len(Con_Const)):
if Con_Const[i] < 0:
Con_Const[i] = -1*Con_Const[i];
Con_Coe[i] = list(-1*vector(Con_Coe[i]));
Constraint_Type[i] = flip_sign(Constraint_Type[i]);
is_flip = true;
if is_flip == true:
P = InteractiveLPProblem(Con_Coe, Con_Const, Obj_Coe, Variables_Names, problem_type=Problem_Type,constraint_type=Constraint_Type,variable_type=Variable_Type,);
show('Slack, Surplus and Artificial Variables are constructed')
#initialize artificial, slack and surplus variables_list
artificial_variables_list = [];
slack_variables_list = [];
surplus_variables_list = [];
slack_surplus_variables_list = [];
#create artificial variables
# initialize counter "l"
l = 1;
#add artificial variables for constraints ">=" and "=="
for i in range(len(Constraint_Type)):
if Constraint_Type[i] == "==" or Constraint_Type[i] == ">=":
artificial_variables_list += [["R_"+str(l),i,1]];
l += 1;
#add slack and surplus variables
l = 1;
for i in range(len(Constraint_Type)):
if Constraint_Type[i] == "<=":
slack_variables_list += [["s_"+str(l),i,1]];
slack_surplus_variables_list += [["s_"+str(l),i,1]];
l += 1;
if Constraint_Type[i] == ">=":
surplus_variables_list += [["s_"+str(l),i,-1]];
slack_surplus_variables_list += [["s_"+str(l),i,-1]];
l += 1;
additional_variables_list = artificial_variables_list + slack_surplus_variables_list;
# make full variables list
full_variables_list = list(variables_list);
for i in range(len(additional_variables_list)):
full_variables_list += [additional_variables_list[i][0]]
#make full constraint coefficients
Full_Con_Mat = [];
for row in Con_Coe:
Full_Con_Mat += [list(row)];
for i in range(len(Constraint_Type)):
for j in range(len(additional_variables_list)):
Full_Con_Mat[i] += [0];
for i in range(len(additional_variables_list)):
Full_Con_Mat[additional_variables_list[i][1]][i+len(variables_list)] = additional_variables_list[i][2]
P0_con_type = ["=="]*len(Constraint_Type);
P0_var_list = list(full_variables_list);
P0_var_type = [">="]*len(full_variables_list);
P0_obj_coe = Obj_Coe + [0]*len(additional_variables_list)
P0 = InteractiveLPProblem(Full_Con_Mat, Con_Const, P0_obj_coe, P0_var_list, problem_type=Problem_Type,constraint_type=P0_con_type,variable_type=P0_var_type,);
show(P0)
show('Phase one')
show(html('Calculate:'))
P1_con_type = ["=="]*len(Constraint_Type);
P1_var_list = list(full_variables_list);
P1_var_type = [">="]*len(full_variables_list);
# construct objective coefficient of phase 1
P1_obj_coe = [0]*len(full_variables_list);
for i in range(len(artificial_variables_list)):
P1_obj_coe[len(variables_list)+i] = 1;
P1_prob_type = "min";
P1 = InteractiveLPProblem(Full_Con_Mat, Con_Const, P1_obj_coe, P1_var_list, problem_type=P1_prob_type,constraint_type=P1_con_type,variable_type=P1_var_type,);
show(P1)
#confirm the full variables list for phase one
P1_variables_list = list(P1.x())
# convert P1 to -max problem
P1_prob_type = "-max"
for i in range(len(artificial_variables_list)):
P1_obj_coe[len(variables_list)+i] = -1;
P1 = InteractiveLPProblem(Full_Con_Mat, Con_Const, P1_obj_coe, P1_var_list, problem_type=P1_prob_type,constraint_type=P1_con_type,variable_type=P1_var_type,);
#show(P1)
# revise the objective function
show("row operation is done to the r-row to resolve inconsistency")
P1_revobj_coe = [0]*len(full_variables_list);
for vr in artificial_variables_list:
P1_revobj_coe = list(vector(P1_revobj_coe) + vector(Full_Con_Mat[vr[1]]));
P1_revobj_coe = list(vector(P1_revobj_coe) + vector(P1_obj_coe))
# calculate the revised objective_constant
P1_revobj_const = 0;
for vr in artificial_variables_list:
P1_revobj_const = P1_revobj_const - Con_Const[vr[1]];
#show min problem
P1 = InteractiveLPProblem(Full_Con_Mat, Con_Const, list(-1*vector(P1_revobj_coe)), P1_var_list, problem_type='min',constraint_type=P1_con_type,variable_type=P1_var_type,objective_constant_term = -1*P1_revobj_const);
show(P1)
# change back to -max
P1 = InteractiveLPProblem(Full_Con_Mat, Con_Const, P1_revobj_coe, P1_var_list, problem_type=P1_prob_type,constraint_type=P1_con_type,variable_type=P1_var_type,objective_constant_term = P1_revobj_const);
#construct dictionary for a packaged problem
#construct coefficient matrix without column of artificial variables and slack variables
trim_Con_Mat = [];
for row in Con_Coe:
trim_Con_Mat += [list(row)];
for i in range(len(Con_Const)):
for j in range(len(surplus_variables_list)):
trim_Con_Mat[i] += [0];
for j in range(len(surplus_variables_list)):
trim_Con_Mat[surplus_variables_list[j][1]][len(variables_list)+j] = -1;
trim_obj_coe = [0]*len(trim_Con_Mat[0]);
for var in artificial_variables_list:
for j in range(len(trim_Con_Mat[0])):
trim_obj_coe[j] += trim_Con_Mat[var[1]][j];
#get the decision and surplus variables list
vr_list = list(variables_list);
for vr in surplus_variables_list:
vr_list += [vr[0]];
#construct artificial and slack variables list according to order of constaints
art_slk_variables_list = [0]*len(artificial_variables_list + slack_variables_list);
for vr in artificial_variables_list:
art_slk_variables_list[vr[1]] = vr[0]
for vr in slack_variables_list:
art_slk_variables_list[vr[1]] = vr[0]
P1 = InteractiveLPProblemStandardForm(trim_Con_Mat, Con_Const, trim_obj_coe, x =vr_list, slack_variables = art_slk_variables_list, objective_constant_term = P1_revobj_const,problem_type = "-max",objective_name = "-r")
D1 = P1.initial_dictionary()
show(D1)
# solve phase one
iter = 1;
while D1.is_optimal()==false:
show(html('Iteration step ' + str(iter) + '.'))
dis_tableau(tableau(D1,P1_variables_list,'min','r'))
OC_list=list((D1. objective_coefficients()));
enter_index = OC_list.index(max(OC_list));
D1.enter(D1.nonbasic_variables()[enter_index]);
show(html("Enter " + str(D1.entering()) ))
show(D1)
show(table(ratio_table(D1)))
entered_ratio=D1.ratios();
possible_leaving_variables = [row[1] for row in entered_ratio];
possible_leaving_ratios = [row[0] for row in entered_ratio];
leave_index = possible_leaving_ratios.index(min(possible_leaving_ratios));
D1.leave(possible_leaving_variables[leave_index]);
show(html("Leave " + str(D1.leaving()) ))
show(D1)
D1.update();
iter += 1;
dis_tableau(tableau(D1,P1_variables_list,'min','r'))
show(D1)
#check feasibility
if D1.objective_value() != 0:
show('Min r \\$\\neq\\$ 0, the LP problem is not feasible.')
assert()
#run phase two if ok
show(html('
Minimum of sum of artificial variables is zero, proceed to Phase two
'))
show('Phase two')
#retrieve constraint constants from D1
P2_con_const = list(D1.constant_terms());
#construct variable lists for phase two
P1f_variables_list = list(P0.x());
P2_variables_list = [];
for i in range(len(variables_list)):
P2_variables_list += [P1f_variables_list[i]];
for i in range(len(slack_surplus_variables_list)):
P2_variables_list += [P1f_variables_list[len(variables_list)+len(artificial_variables_list)+i]];
#retrieve variable order from D1
D1_basic_variables_list = list(D1.basic_variables());
D1_nonbasic_variables_list = list(D1.nonbasic_variables());
P2_Con_Mat = [[0]*len(P2_variables_list) for i in range(len(D1_basic_variables_list))];
i = 0;
j = 0;
for vri in D1_basic_variables_list:
j = 0;
for vrj in P2_variables_list:
if vri in D1_basic_variables_list and (vri == vrj):
P2_Con_Mat[i][j] = 1;
if vrj in D1_nonbasic_variables_list:
P2_Con_Mat[i][j] = D1.column_coefficients(vrj)[i]
j += 1;
i += 1;
# construct phase 2 objective coefficient
P2_obj_coe = list(Obj_Coe) + [0]*len(slack_surplus_variables_list)
# set con type
P2_con_type = ["=="]*len(P2_con_const);
P2_var_type = [">="]*len(P2_variables_list);
P2 = InteractiveLPProblem(P2_Con_Mat, P2_con_const, P2_obj_coe, P2_variables_list, problem_type=Problem_Type,constraint_type=P2_con_type,variable_type=P2_var_type,);
if Problem_Type =='-max':
P22_obj_coe = list(-1 * vector(P2_obj_coe));
P22 = InteractiveLPProblem(P2_Con_Mat, P2_con_const, P22_obj_coe, P2_variables_list, problem_type='min',constraint_type=P2_con_type,variable_type=P2_var_type,);
show(P22)
else:
show(P2)
show("row operation is done to the objective function to resolve inconsistency")
# do row operation to objective function
P2_obj_const = 0
i = 0;
for vr in D1_basic_variables_list:
Multiplier = P2_obj_coe[P2_variables_list.index(vr)];
P2_obj_const = P2_obj_const + Multiplier*P2_con_const[i]
for j in range(len(P2_variables_list)):
P2_obj_coe[j] = P2_obj_coe[j] - Multiplier*P2_Con_Mat[i][j]
i += 1;
P2 = InteractiveLPProblem(P2_Con_Mat, P2_con_const, P2_obj_coe, P2_variables_list, problem_type=Problem_Type,constraint_type=P2_con_type,variable_type=P2_var_type,objective_constant_term = P2_obj_const);
if Problem_Type =='-max':
P22_obj_coe = list(-1 * vector(P2_obj_coe));
P22_obj_const = -1*P2_obj_const;
P22 = InteractiveLPProblem(P2_Con_Mat, P2_con_const, P22_obj_coe, P2_variables_list, problem_type='min',constraint_type=P2_con_type,variable_type=P2_var_type,objective_constant_term = P22_obj_const);
show(P22)
else:
show(P2)
# construct dictionary for P2
# construct list of basic and nonbasic variables
P2_init_base_var = list(D1_basic_variables_list);
P2_init_nonbase_var = [];
for vr in D1_nonbasic_variables_list:
if vr in P2_variables_list:
P2_init_nonbase_var +=[vr];
#construct dictionary matrix from D1
D2_init_dict_Mat = list([list([]) for i in range(len(P2_init_base_var))]);
for i in range(len(P2_init_base_var)):
for vr in P2_init_nonbase_var:
D2_init_dict_Mat[i] += [D1.column_coefficients(vr)[i]];
#construct objective coefficients of non basic variables
D2_obj_coe = [];
for vr in P2_init_nonbase_var:
D2_obj_coe += [P2_obj_coe[P2_variables_list.index(vr)]];
#construct variable string
D2_var_string = "";
for i in range(len(P2_init_nonbase_var)):
if i ==0:
D2_var_string += str(P2_init_nonbase_var[i]);
else:
D2_var_string = D2_var_string + "," + str(P2_init_nonbase_var[i]);
for i in range(len(P2_init_base_var)):
D2_var_string = D2_var_string + "," + str(P2_init_base_var[i]);
#construct D2
D2_obj_name = "z"
if Problem_Type =="-max":
D2_obj_name = "-z"
A = matrix(Base_Ring, D2_init_dict_Mat);
b = vector(Base_Ring, P2_con_const);
c = vector(Base_Ring, D2_obj_coe);
R = PolynomialRing(Base_Ring, D2_var_string, order="neglex");
from sage.numerical.interactive_simplex_method import LPDictionary
D2 = LPDictionary(A, b, c, P2_obj_const, R.gens()[len(P2_init_nonbase_var):], R.gens()[:len(P2_init_nonbase_var)],D2_obj_name)
show(D2)
iter = 1
while D2.is_optimal()==false:
show(html('Iteration step ' + str(iter) + '.'))
dis_tableau(tableau(D2,P2_variables_list,Problem_Type,'z'))
OC_list=list((D2. objective_coefficients()));
enter_index = OC_list.index(max(OC_list));
D2.enter(D2.nonbasic_variables()[enter_index]);
show(html("Enter " + str(D2.entering()) ))
show(D2)
show(table(ratio_table(D2)))
entered_ratio=D2.ratios();
possible_leaving_variables = [row[1] for row in entered_ratio];
possible_leaving_ratios = [row[0] for row in entered_ratio];
leave_index = possible_leaving_ratios.index(min(possible_leaving_ratios));
D2.leave(possible_leaving_variables[leave_index]);
show(html("Leave " + str(D2.leaving()) ))
show(D2)
D2.update();
iter += 1;
dis_tableau(tableau(D2,P2_variables_list,Problem_Type,'z'))
show(D2)
#retrieve solution
D2_final_basic_var = list(D2.basic_variables());
D2_final_nonbasic_var = list(D2.nonbasic_variables());
sol = [];
for vr in variables_list:
if vr in D2_final_basic_var:
sol += [D2.constant_terms()[D2_final_basic_var.index(vr)]];
if vr in D2_final_nonbasic_var:
sol += [0];
opt_val = D2.objective_value();
if Problem_Type =="-max":
opt_val = -opt_val;
show("The optimal value is " + str(opt_val) + ", with an optimal solution " + str(vector(sol)))