-1

Please tell me how can I solve my problem with the below variables.

Which solver is better?

Objective function = x^T Q x + C^T x

  • Q is diagonal matrix with Q(i,i)>=0
  • x.shape=(n,1)

Equal constraint : Aeq.x =beq

  • Aeq.shape=(n,n) and beq.shape=(n,1)

Lower bound and upper bound : lb <= x <= ub

import random
import numpy as np

lb.shape=(n,1) and ub.shape=(n,1)

I=np.eye(24)
Z=np.zeros((24,24))
a=0.012
b=1.1
gamma1=0.9/80
gamma2=1.1/80
MM=np.eye(24)

for i in range (22):
    MM[i+1,i]=-1
    MM[0,23]=-1

M=random.randint(200,300, size=(24,1))
max_pch=30.0
max_pdch=30.0
ppp=random.randint(150,200, size=(24,))
Q=np.asarray(np.bmat([[a*I,Z,Z,Z],[Z,a*I,Z,Z],[Z,Z,0.00001*I,Z],[Z,Z,Z,0.00001*I]  ]))
C=np.asarray(np.bmat([[b*np.ones(24),b*np.ones(24),0*np.ones(24),ppp]]))

Aeq=np.asarray(np.bmat([[-I,I,Z,I], [-gamma1*I, gamma2*I,MM,Z],[Z,Z,Z,Z],[Z,Z,Z,Z]]))
beq=np.asarray(np.bmat([[M],[np.zeros((72,1))]]))

lb=np.asarray(np.bmat([[0*np.ones(24),0*np.ones(24),[0.1],0.1*np.ones(22),[0.9],0.0*np.ones(24)]]))

ub=np.asarray(np.bmat([[max_pch*np.ones(24),max_pdch*np.ones(24),[0.9],0.9*np.ones(22),[0.9],500*np.ones(24)]]))

x = solve_qp(P=Q, q=C.T.reshape((96,)),
                G=None , h=None,
                A=Aeq , b=beq.reshape(96,),
                lb=lb.T.reshape((96,)) , ub=ub.T.reshape((96,)))
print("QP solution: x = {}".format(x))

What is the problem?

  • QP solution: x = None

The same code in Matlab (with fmincon) gives the correct result. However, in Python, I can't get that result.

Alex Waygood
  • 6,304
  • 3
  • 24
  • 46
sajad parvizi
  • 1
  • 2

2 Answers2

0

If you print Q you will see it has rows and columns that are all zeros. This makes the matrix Q no longer strictly positive definite (PD). You can either use a solver that allows for positive semi-definite (PSD) matrices (most QP solvers will allow this, just not the one you are using now), or you can add a small number (say 1e-6) to all zero diagonal entries.

Erwin Kalvelagen
  • 15,677
  • 2
  • 14
  • 39
0

It seems you are using the solve_qp function from qpsolvers. From the FAQ of the package's README:

I have a non-convex quadratic program. Is there a solver I can use?

  • Unfortunately most available QP solvers are designed for convex problems.
  • If your cost matrix P is semi-definite rather than definite, try OSQP.

I tried your problem with OSQP, and the solver exited with the status "primal infeasible". This means the solver was able to solve your problem and found a certificate that it has no solution. (I don't know what fmincon returned, but you will want to check that its results satisfies all constraints Aeq * x == beq and lb <= x <= ub.)

Other advice / things you can try:

  • The last lines of Aeq and beq are zeros. In that case, it's best not to add these lines to the problem. Some solvers don't handle 0 * x == 0 lines well.
  • Play with enabling/disable constraints to get a feeling of your problem. For example, if we disable lb, we see that the values in the middle of the solution vector hit 30 and that forces the other values to become <= -500. Either the objective is steering there or this is caused by Aeq * x == beq (most likely the latter since the problem with lb seems infeasible).
  • You can relax the hard equality constraint by converting it to a quadratic objective weight * || Aeq * x - beq ||^2 in the objective function. This way the problem will always be feasible, and you can play with the weight parameter to understand the effect of this constraint.

Good luck with your study!

Tastalian
  • 369
  • 4
  • 12