Simple Inequality Constraint Breaks QP Problem (Control Allocation)

Question

Optimization problem: Just for context, not strictly needed. Using SciPy's "minimize" the code is allocating optimal control forces to thrusters on a ship given environmental forces acting on the ship (wind, waves & current) in the 3 DOFs surge, sway and yaw to hold the ship's position (Important for Dynamic Positioning). This ship has 4 thrusters, 1 fixed-angle thruster (tunnel) and 3 rotating thrusters (azimuth). The fixed thruster has a force in one direction (x1) and the rotating thrusters has forces in two directions, x and y, implemented as variable pairs (x2, x3), (x4, x5), (x6, x7). E.g. x2 and x3 is the force in x and y, respectively, for rotating thruster 1. There are also 3 slack variables for the 3 DOFs to ensure a solution.

Problem on standard form

Programming problem: I want to limit one of the thrusters' force by implementing an inequality constraint on it, but by implementing the constraint "ineq1" all outputs are minimized to zero. The constraint should limit x1 and is very simply like this: (x1 < x1_max) Follow-up: "minimize" seems to ignore the first constraint in the constraint dictionary and it doesn't matter if it's an equality or inequality constraint. The second constraint is however upheld.

Output with only equality constraint: Correct output

Output with breaking inequality constraint: Incorrect output with inequality

What I've tried is

• ... trying different algorithms
• ... test every combination of positives/negatives in the constraint
• ... use an array/vector A as in the standard form A*x < b instead of just x[0]
• ... using a number instead of fmax[0]
• ... tried other design variables and combinations of them

Code:

This can be readily run. Change between "cons1" and "cons2" in line 42 to see the different results.

import numpy as np
from scipy.optimize import minimize

# Parameters
nf  = 1         # No. of fixed thrusters
nr  = 3         # No. of rotating thrusters
nt  = nf + 2*nr # No. of actuators
ns  = 3         # No. of slack variables (slacks)
n   = nt + ns   # Design variables
Lx  = np.array([30.5, 22.25, -35.5, -35.5])         # Thruster position
W   = np.diag([1, 1, 1, 1, 1, 1, 1])                # Weighting for thrusters
Q   = np.diag([1, 1, 1])*5000                       # Weighting for slacks
Te  = np.array([[0,1,0,1,0,1,0],                    # Thrust config. matrix
                [1,0,1,0,1,0,1],
                [Lx[0],0,Lx[1],0,Lx[2],0,Lx[3]]])
fmax = np.array([ 50000,  60000,  75000,  75000])

tau = np.array([50000,35000,0]) # Environmental forces in surge, sway & yaw
x0 = np.zeros(n)

# Problem definition
def obj(x):
    fe = x[:nt]
    s  = x[nt:]
    return fe.transpose()@W@fe + s.transpose()@Q@s

def eqcon(x):
    fe = x[:nt]
    s  = x[nt:]
    return tau + s - Te@fe # = 0

def ineq1(x):
    return fmax[0] - x[0] # > 0

cons1 = {'type': 'eq',   'fun': eqcon}

cons2 = {'type': 'eq',   'fun': eqcon,
         'type': 'ineq', 'fun': ineq1}

# Computing
opt = {'disp': True}
res = minimize(obj, x0, method='SLSQP', constraints=cons2, options=opt)

# Debugging
print("Environmental forces (surge, sway & yaw): ", tau)
print("Forces: ", np.round(res.x[:nt],1))
print("Slacks: ", np.round(res.x[nt:],1))
print(f"Thruster 1 has angle 90.00 with force {res.x[0]/1000:0.2f} kN")

angles = np.zeros(nr)
ux = res.x[nf:nt][::2]      # Every other force element in x for azimuth: ux2, ux3, ...
uy = res.x[nf:nt][1:][::2]  # uy2, uy3, ...

for i in range(nr):
    angles[i] = np.arctan2(uy[i], ux[i]) * 180 / np.pi
    print(f"Thruster {nf+i+1} has angle {angles[i]:0.2f} with force {np.sqrt(ux[i]**2+uy[i]**2)/1000:0.2f} kN")

Any ideas?

Answer 1

The constraints are expected to be a list of dictionaries, so using

cons2 = [{'type': 'eq',   'fun': eqcon},
         {'type': 'ineq', 'fun': ineq1}]

should fix your problem. That being said, the keys inside a dictionary should be unique. Otherwise, you aren't able to access a specific value in the dictionary by ambiguous keys.