With the help of nbruin, this is a perfect way to resolve a double sum problem.
c = matrix(3,3,(-1,2,-3,3,2,-1,2,1,3))
d = matrix(3,3,(1,2,3,3,2,1,2,1,3))
show(c)
show(d)
p = MixedIntegerLinearProgram(solver='GLPK')
v = p.new_variable(nonnegative=True)
for j in range(3):
p.set_objective(sum(c[i,j]*v[i,j] for i in range(0,3) for j in range(0,3)))
p.add_constraint(sum(d[i]*v[i,j] for i in range(0,3)) <= 1)
p.show()
p.solve()
vv=p.get_values(v)
show(vv)
w=v.items()
w
But 1) w is not sorted 2) I would like to construct a table with $x_{(a, b)}$ in one column and the optimal value of $x_c$ which corresponds to the $x_{(a, b)}$ .