FDA   = 0.0306
rho   = 1.23
mu    = 1.78e-5
SwetS = 2.05
k     = 1.2
e     = 0.96
W0    = 4940
Nlift = 2.5
tau   = 0.12
Vmin  = 22
CLmax = 2.0

import numpy as np

# minimize     c^T x
# subject to   Ax == b
#              Cx <= d
#              log_sum_exp(Ex+f) <= 0
#              log_sum_exp(Gx+h) <= 0
#              log_sum_exp(Px+q) <= 0

import cvxpy as cp

# x: [  a, c_d, c_f, c_l,   r,   s,   v,   w, w_w]
x = cp.Variable(9)

# min c^T x  
             #  a, c_d, c_f, c_l,   r,   s,   v,   w, w_w#
c = np.array([  0,   1,   0,   0,   0,   1,   2,   0,   0])    # s +c_d +2v
obj = cp.Minimize(c.T@x)

# Ax == b 
              #  a, c_d, c_f, c_l,   r,   s,   v,   w, w_w#
A = np.array([[  0,   0,   1,   0,  .2,   0,   0,   0,   0],  # c_f +.2r
              [  0,   0,   0,  -1,   0,  -1,  -2,   1,   0],  # w -2v -c_l -s
              [-.5,   0,   0,   0,  -1,  .5,   1,   0,   0]]) # v -r +.5s -.5a
b = -np.array([-np.log(0.074), np.log(2/rho), np.log(rho/mu)])

# Cx <= d
             #  a, c_d, c_f, c_l,   r,   s,   v,   w, w_w#
C = np.array([  0,   0,   0,   0,   0,  -1,   0,   1,   0])   # w -s
d = -np.log(2/rho/Vmin/Vmin/CLmax)

# log_sum_exp(Ex+f) <= 0
              #  a, c_d, c_f, c_l,   r,   s,   v,   w, w_w#
E = np.array([[  0,  -1,   0,   0,   0,  -1,   0,   0,   0],   # -c_d -s
              [  0,  -1,   1,   0,   0,   0,   0,   0,   0],   # c_f -c_d
              [ -1,  -1,   0,   2,   0,   0,   0,   0,   0]])  # 2c_l -c_d -a
f =  np.array([np.log(FDA), np.log(k*SwetS), -np.log(np.pi*e)])

# log_sum_exp(Gx+h) <= 0
              #  a, c_d, c_f, c_l,   r,   s,   v,   w, w_w#
G = np.array([[  0,   0,   0,   0,   0,   0,   0,  -1,   0],   # -w
              [  0,   0,   0,   0,   0,   0,   0,  -1,   1]])  # w_w -w
h =  np.array([np.log(W0), 0])

# log_sum_exp(Px+q) <= 0
              #  a, c_d, c_f, c_l,   r,   s,   v,   w, w_w#
P = np.array([[  0,   0,   0,   0,   0,   1,   0,   0,  -1],  # s -w_w
              [1.5,   0,   0,   0,   0,  .5,   0,  .5,  -1]]) # 1.5a + .5w + .5s -w_w
q =  np.array([np.log(45.42), np.log(8.71e-5*Nlift*np.sqrt(W0)/tau)])

constr = [A@x == b, C@x <= d, 
          cp.log_sum_exp(E@x+f) <= 0, 
          cp.log_sum_exp(G@x+h) <= 0, 
          cp.log_sum_exp(P@x+q) <= 0]

cp.Problem(obj, constr).solve(verbose=True)

AR  = np.exp(x.value[0])
C_D = np.exp(x.value[1])
C_f = np.exp(x.value[2])
C_L = np.exp(x.value[3])
Re  = np.exp(x.value[4])
S   = np.exp(x.value[5])
V   = np.exp(x.value[6])
W   = np.exp(x.value[7])
W_w = np.exp(x.value[8])

print(f'AR  = {AR}')
print(f'S   = {S}')
print(f'C_D = {C_D}')
print(f'C_L = {C_L}')
print(f'C_f = {C_f}')
print(f'Re  = {Re}')
print(f'W   = {W}')
print(f'W_w = {W_w}')
print(f'V   = {V}')