CVXOPT with only equality constraints

Question

I am trying the following as learning exercise in CVXOPT. I have made minor modifications to the example code here by removing the inequality constraints and adding few more equality constraints.

from cvxopt import solvers, blas, matrix, spmatrix, spdiag, log, div
solvers.options['show_progress'] = False
import numpy as np    
np.random.seed(1)

# minimize     p'*log(p)
# subject to
#              sum(p) = 1
#              sum(p'*a) = target1
#              sum(p'*max(a-K,a^2)) = target2

a = np.random.randint(20, 30, size=500)
target1 = 30
target2 = 0.60
K = 26

A = matrix(np.vstack([np.ones(500), a, np.array([max(x-K,x*x) for x in a])]))
b = matrix([1.0, target1, target2])

n = 500
def F(x=None, z=None):
   if x is None: return 0, matrix(1.0, (n,1))
   if min(x) <= 0: return None
   f = x.T*log(x)
   grad = 1.0 + log(x)
   if z is None: return f, grad.T
   H = spdiag(z[0] * x**-1)
   return f, grad.T, H
sol = solvers.cp(F, A=A, b=b)
p = sol['x']

But when I perform the following:

np.sum(p)
243.52686763225338

This violated the first constraint of the optimization. I am not able to figure what is going wrong here. (Please note since I am using random numbers to generate variable a your np.sum(p) will produce different values but you should observe the same violation as mine.

Even if I keep the inequality constraints from the original link and add the two additional equality constraints, the equality constraints are violated.

Is there any other package I can use reliably i.e a package which is maintained?

Edit: If there is no feasible solution, shouldn't there be a message that no feasible solution found?

Answer 1

As @tihom commented the problem is infeasible. Are you really sure this is the problem you want to solve? Your first constraint implies:

p1 + p2 + ... + pn = 1
p1*a1 + p2*a2 + ... + an*pn = 30
p1*a1^2 + p2*a2^2 + ... pn*an^2 = 0.6

The last constraint can never be satisfied at the same time as the first or the second one since ai >= 20 for all i. That is, the sum p1*a1^2 + p2*a2^2 + ... pn*an^2 is always greater than the other sums (note that pi > 0).

If you let target1 = sum(a/500.) and target2= sum(a*a/500.) there exists a point satisfying your constraints and you can find the optimal solution.

Note that in the last constraint the maximum simplifies to max(a - K, a^2) = a^2, which is true regardless of a.

Edit: If you inspect the solution (e.g., print sol) you will get something like:

{'status': 'unknown', 'zl': <0x1 matrix, tc='d'>, 'dual slack': 1.0000000000000007, 'relative gap': 0.005911420508296136, 'dual objective': -97.17320604198335, 'snl': <0x1 matrix, tc='d'>, 'gap': 0.9737154924375709, 'primal objective': -164.7176835197311, 'primal slack': 0.9737154924375703, 'znl': <0x1 matrix, tc='d'>, 'primal infeasibility': 0.5114570271204905, 'dual infeasibility': 0.5091221046374248, 'sl': <0x1 matrix, tc='d'>, 'y': <3x1 matrix, tc='d'>, 'x': <500x1 matrix, tc='d'>}

Notice that the status is 'unknown', i.e., no feasible solution is found. This is in the documentation of cvxopt.solvers.cp: http://cvxopt.org/userguide/solvers.html?highlight=cp#cvxopt.solvers.cp