0

I am trying to solve an optimization problem using my own code. I want to use Sequential Quadratic Programming to minimize a n objective function given a few inequality constraints. I have made an active-set solver for the sub-QP problem and I have tested it out. The problem is when I have a non-quadratic function and I try to make approximizations at each iteration so that I can solve multiple sub-QP problems to find the solution, my first iteration (Which I'm doing by hand right now) is returning a point that is feasible in the linearized approximations, however it is not feasible for the original non-linear constraints. This means that my algorithm cannoot proceed as the sub-QP problem requires a feasible point!

This is my code

x0= np.array([[0.1], [0.2]])
x = x0.copy()
H = np.identity(len(x0))
cons = [g1,g2,g3,g4,g5,g6]

# Hessian Initial guess 
H = np.identity(len(x0))
c = grad(f,x0).T - 0.5*(x0.T@H.T) - 0.5*(x0.T@H)
A = np.zeros([len(cons),len(x0)])
b = np.zeros((len(cons),1))
it=0
for i in range(len(cons)):
  A[i,:] = grad(cons[i],x0).T
  b[i] = (cons[i](x0) - grad(cons[i],x0).T@x0)

x = qpcon(H,c.T,A,b,x0,1e-6)
print(g1(x))

The qpcon function is my own active set code:

# Equality Constrained (Needed for the inequality constrained)
def qp(H,c,A,b):
    Qinv = H
    lam = -np.linalg.inv(A@Qinv@A.T)@A@Qinv@c + np.linalg.inv(A@Qinv@A.T)@b
    x = -Qinv@A.T@lam - Qinv@c
    return x, lam


#Inequality Constrained
def qpcon(H,q,C,d,x0,epsilon):
  Q = np.linalg.inv(H)
  k = 0
  xk = x0
  w_check = abs(A@x0 + b) <1e-6*np.ones([len(b),1])
  Wk = np.array([],dtype=int)
  for i in range(len(w_check)):
    if w_check[i] == 1:
      Wk = np.append(Wk,i)
  while True:
    Cw = np.zeros([len(Wk),len(x0)])
    dw = np.zeros([len(Wk),1])
    Cw = C[Wk]
    dw = d[Wk]
    if len(Wk) == 0:
      p =qpu(H,q+Q.T@xk)
      sigma = 0
    else:
      p, sigma = qp(H,q+Q.T@xk,Cw,np.zeros((len(Wk),1)))
    if np.linalg.norm(p)<epsilon:
      if (sigma>=0).all():
        xstar = xk
        break
      else:
        i = np.argmin(sigma)
        Wkp1 = np.delete(Wk,i)
        xkp1 = xk
    else:
      alpha = 1
      B= np.array([])
      for i in np.setdiff1d(np.arange(len(w_check)), Wk, assume_unique=False):
        if C[i].T@p>0:
          alpha_b = -(C[i].T@xk+d[i])/(C[i].T@p)
          if alpha_b<alpha:
            alpha = alpha_b
            B = i
      Wkp1 = np.append(Wk,B)
      Wkp1 = np.array([ int(x) for x in Wkp1])
      xkp1 = xk + alpha*p
    Wk = Wkp1.copy()
    xk = xkp1.copy()
    k = k + 1
  return xstar

The Objective function and constrains are as follows (It's a beam problem where the volume has to be minimized but it is also constrained by the maximum Von-Mises stress and maximum deflection, g3(x) through g4(x) are just bounds on the design variables. All non-equalities are in the form g(x)<= 0 )

def f(x):
    L = 1
    w = x[0]
    h = x[1]
    return L*w*h

def g1(x):
    w = x[0]
    h = x[1]
    L = 1
    S_y = 400e6
    n = 2
    F = 1e5
    return (6*F*L/(w*h**2) - S_y/n)

def g2(x):
    w = x[0]
    h = x[1]
    E = 210e9
    L = 1
    F = 1e5
    d = 1e-2
    return ((4*F*L**3)/(E*w*h**3) - d)

def g3(x):
    w = x[0]
    h = x[1]
    return x[0] - 400e-3

def g4(x):
    w = x[0]
    h = x[1]
    return x[1] - 400e-3

def g5(x):
    w = x[0]
    h = x[1]
    return 0.2e-3 - x[0]

def g6(x):
    w = x[0]
    h = x[1]
    return 0.2e-3 - x[1]

Please let me know if I'm doing something wrong or if there is a simple way to remedy this infeasible step.

I have tried the "HANDLING INCONSISTENT LINEARIZATIONS" method from section 18.3 of book Numerical Optimization by Jorge Nocedal & Stephen J. Wright.

Mohammad
  • 1

0 Answers0