Optimizing with SciPy

Before we can use SciPy, we have to install it. Thankfully, SciPy is one of the best supported python packages in the world, so to install it one can just use pip install scipy.

SciPy has many optimizers available, but I'm going to focus on using SLSQP since it can use the gradient information from OptVL and supports constraints. To use SciPy's SLSQP, we will need to supply it with custom objective and constraint functions as well as the derivatives for each. These functions need to take in the design variables and apply them to our OptVL solver. The snippet below provides an example of an objective function.

 def objective_function(x):
    ovl.set_constraint("Elevator", x[0])
    ovl.set_surface_params({"Wing":{"aincs":x[1:]}})

    ovl.execute_run()
    cd = ovl.get_total_forces()['CD']
    print(x, cd)

    return cd

Note: the objective function is specified by the return value of the function. If you wanted to save data about each iteration or write output, the objective function would be a good place to add that functionality. To supply the gradient information to SLSQP, we have to define another function that returns the gradients for a given design variable vector.

def objective_gradient(x):
    # Partial derivatives of the objective_function
    ovl.set_constraint("Elevator", x[0])
    ovl.set_surface_params({"Wing":{"aincs":x[1:]}})

    ovl.execute_run()

    sens = ovl.execute_run_sensitivities(['CD'])
    dcd_dele = sens['CD']['Elevator']
    dcd_daincs = sens['CD']['Wing']['aincs']

    # concatinate the two and return the derivs
    return np.concatenate(([dcd_dele], dcd_daincs))

The function ovl.execute_run_sensitivities(['CD']) does all the necessary work to compute the derivatives for the given list of functions. We just need to parse the sens dictionary for the derivatives with respect to the design variables we are interested in.

One also needs to repeat the process for the constraints. This involves creating functions for both the constraint values and their corresponding gradients.

Example script

The script below shows the full optimization script for minimizing the drag of the example aircraft in trim, with respect to the wing twist and elevator position.

"""A scipy based optimization to trim an aircraft for an aicraft using optvl"""
import numpy as np
from scipy.optimize import minimize
from optvl import OVLSolver

ovl_solver = OVLSolver(geo_file="aircraft.avl", debug=False)

# setup OptVL 
ovl_solver.set_parameter("Mach", 0.0)

# Define your custom objective function with outputs from OptVL
def objective_function(x):
    ovl_solver.set_constraint("Elevator", x[0])
    ovl_solver.set_surface_params({"Wing":{"aincs":x[1:]}})

    ovl_solver.execute_run()
    cd = ovl_solver.get_total_forces()['CD']
    print(x, cd)

    return cd

def objective_gradient(x):
    # Partial derivatives of the objective_function

    # we are trusting that the design variables have already been applied 
    # and propogated through by the objective_function. 

    ovl_solver.set_constraint("Elevator", x[0])
    ovl_solver.set_surface_params({"Wing":{"aincs":x[1:]}})

    ovl_solver.execute_run()

    sens = ovl_solver.execute_run_sensitivities(['CD'])
    dcd_dele = sens['CD']['Elevator']
    dcd_daincs = sens['CD']['Wing']['aincs']

    # concatinate the two and return the derivs
    return np.concatenate(([dcd_dele], dcd_daincs))


# Define equality constraint: h(x) = 0
def eq_constraint(x):
    ovl_solver.set_constraint("Elevator", x[0])
    ovl_solver.set_surface_params({"Wing":{"aincs":x[1:]}})

    ovl_solver.execute_run()

    # the objective must always be run first
    coeff = ovl_solver.get_total_forces()

    cl_target = 1.5
    cm_target = 0.0

    cl_con = coeff['CL'] - cl_target
    cm_con = coeff['CM'] - cm_target

    return np.array([cl_con, cm_con])

# Define the gradient of the equality constraint
def eq_constraint_jac(x):
    sens = ovl_solver.execute_run_sensitivities(['CL', 'CM'])
    dcl_dele = sens['CL']['Elevator']
    dcl_daincs = sens['CL']['Wing']['aincs']
    dcl_dx = np.concatenate(([dcl_dele], dcl_daincs))

    dcm_dele = sens['CM']['Elevator']
    dcm_daincs = sens['CM']['Wing']['aincs']
    dcm_dx = np.concatenate(([dcm_dele], dcm_daincs))

    # concatinate the two and return the derivs
    return np.array([dcl_dx, dcm_dx])

num_sec = 5
# Initial guess for the variables
x0 = np.zeros(1+num_sec)

# Define constraints with their gradients
constraints = [
    {'type': 'eq', 'fun': eq_constraint, 'jac': eq_constraint_jac},   # Equality constraint
]

# Call the minimize function with constraints and their gradients
result = minimize(
    fun=objective_function,       # The objective function to minimize
    jac=objective_gradient,       # The gradient function of the objective
    x0=x0,                     # Initial guess
    constraints=constraints,   # Constraints with gradients
    method='SLSQP',            # Optimization method that supports constraints
    options={'disp': True},     # Display convergence messages
    tol=1e-10
)

# Print the result
print("Optimization result:")
opt_x = result.x
print(f"Elevator: {opt_x[0]}")
print(f"Wing twist: {opt_x[1:]}")
print(f"Final CD: {result.fun}")          # Final value of the objective function

Which should produce the output after the optimization

Optimization terminated successfully    (Exit mode 0)
            Current function value: 0.08551861616299436
            Iterations: 56
            Function evaluations: 58
            Gradient evaluations: 56
Optimization result:
Elevator: 0.09925034689068805
Wing twist: [11.51424344  2.00085521  2.20301866  2.32160102 -4.79748301]
Final CD: 0.08551861616299436