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?
As @tihom commented the problem is infeasible. Are you really sure this is the problem you want to solve? Your first constraint implies:
The last constraint can never be satisfied at the same time as the first or the second one since
ai >= 20
for alli
. That is, the sump1*a1^2 + p2*a2^2 + ... pn*an^2
is always greater than the other sums (note thatpi > 0
).If you let
target1 = sum(a/500.)
andtarget2= 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 ofa
.Edit: If you inspect the solution (e.g.,
print sol
) you will get something like:Notice that the
status
is'unknown'
, i.e., no feasible solution is found. This is in the documentation ofcvxopt.solvers.cp
: http://cvxopt.org/userguide/solvers.html?highlight=cp#cvxopt.solvers.cp