Optimization with python

5.2. Optimization with python#

Written by Marc Budinger (INSA Toulouse), Toulouse, France

The objective of this notebook is to learn how to implement a sizing code and use a simple numerical optimization into a python notebook. The system studied a Electromechanical Actuator (EMA) used for Thrust Vector Control (TVC) of a launcher.

For interested readers, more information can be found in :

the dedicated laboratories of this Jupyter book
in the following document

Delbecq, S., Budinger, M., & Reysset, A. (2020). Benchmarking of monolithic MDO formulations and derivative computation techniques using OpenMDAO. Structural and Multidisciplinary Optimization, 62(2), 645-666. Link

5.2.1. Objectives and specifications#

We want to optimize the selection of a TVC motor/reducer set: dynamics are high and motor inertia effect is not negligeable.

EMA

The sizing scenario is a sinusoïdal displacement with a constant force and an inertia effect:

Specification
static force	40 kN
Inertia	800 kgg
Displacement magnitude	10 mm
Frequency	5 Hz

We first define the specifications and assumptions for the sizing:

from math import pi

# Specifications
Xmax = 10e-3  # [m]
f = 5  # [Hz]
w = 2*pi*5 # [rad/s] angular frequency

Xpmax = w*Xmax  # [m/s] max linear speed
Xppmax = w**2*Xmax  # [m/s²] max linear acceleration

M = 800 # [kg] Equivalent mass inertia of the nozzle
Fstat = 40e3 # [N] static force

A reference example of motor is:

Characteristics
Nom. torque	3.14 N.m
Max. torque	13.4 N.m
Max. speed	7200 rpm
Inertia	2.90E-04 kg.m²
Mass	3.8 kg

We assume here that :

the pitch \(p\) of the roller screw is 10 mm/rev.
The efficiency of mechanical power transmission is 100%.

The objective is to select the reduction ratio of the gear reducer in order to minimize the mass of the motor.

We then define the main characteristics of components and the reference parameters for the scaling laws:

# Assumptions
pitch = 10e-3 / 2 / 3.14  # [m/rad]
nu_screw = 1  # [-]

# Reference motor
T_mot_nom_ref = 3.14 # [N.m] Max torque
T_mot_max_ref = 13.4  # [N.m] Rated torque
W_mot_max_ref = 754  # [rad/s] max speed
J_mot_ref = 2.9e-4 / 2  # [kg.m²] Inertie 
M_mot_ref = 3.8  # [kg] Mass

5.2.2. Objective and constraints#

The objective is to select the motor and the reduction ratio of a gear reducer in order to minimize the mass of the motor.

The design variables are :

the reduction ratio of the reducer
the selection torque of the motor

The objective is the motor mass.

The constraints which should be positives are here:

the speed margin, ie. the motor doesn’t exceed its maximum speed
the torque margin, ie. the motor doesn’t exceed its maximum torque

5.2.3. Sizing code#

The sizing code is defined here in a function which can give an evaluation of the objective and of the constraints function of design variables.

Exercise 5.2

Give equations which link forces and speeds at nozzle and motor levels.
Give scaling laws usefull for the problem.
Express contraints of the design problem

Solution to

For the sizing scenario :

The maximum load force (nozzle) is : \(F_{max} = F_{stat} + M a_{max}\) with \(a_{max} = max(\ddot{x})\)
The maximum torque at motor level is : \(T_{mot, max} = F_{max} \frac{p}{N} + J_{mot} a_{max} \frac{N}{p}\)

The scaling laws are:

Motor mass : \(M_{mot} = M_{mot,ref} \bigl (\frac{T_{mot}}{T_{mot,nom,ref}} \bigr) ^\frac{3}{3.5}\)
Motor inertia : \(J_{mot} = J_{mot,ref} \bigl (\frac{T_{mot}}{T_{mot,nom,ref}} \bigr) ^\frac{5}{3.5}\)
Max motor speed : \(\omega_{mot,max} = \omega_{mot,max,ref} \bigl (\frac{T_{mot}}{T_{mot,nom,ref}} \bigr) ^\frac{-1}{3.5}\)
Max motor torque : \(T_{mot,max} = T_{mot,max,ref} \bigl (\frac{T_{mot}}{T_{mot,nom,ref}} \bigr) ^1\)

We have to check that:

Max speed possible of the motor (coming from scaling laws) > Max speed of the sizing scenario
\(\omega_{mot,max} > max(\dot{x})\frac{N}{p}\)
Max torque possible of the motor (coming from scaling laws) > Max torque of the sizing scenario
\(T_{mot,max,scaling} > T_{mot,max,scenario}\)

def sizing_code(param, arg="print"):
    #  Design variables
    N = param[0]  # Reduction ratio of the reducer
    T_mot = param[1]  # Motor nominal torque

    # --------------------------------------
    # Parameter estimation with scaling laws

    # Motor
    M_mot = M_mot_ref * (T_mot / T_mot_nom_ref) ** (3 / 3.5)
    J_mot = J_mot_ref * (T_mot / T_mot_nom_ref) ** (5 / 3.5)
    W_mot = W_mot_max_ref * (T_mot / T_mot_nom_ref) ** (-1 / 3.5)
    T_mot_max = T_mot_max_ref * (T_mot / T_mot_nom_ref) ** (1)

    # --------------------------------------
    # Load force
    max_load = Fstat + M*Xppmax   # [N]
    
    # --------------------------------------
    # Motor torque calculation with motor inertia
    T_mot_scenario = (
        max_load * pitch / N / nu_screw
        + J_mot * Xppmax *  N / pitch    ) # [N.m] max electromagnetic torque
    
    # --------------------------------------
    # Objectives and constrants calculations

    # Objective = motor mass
    objective = M_mot

    # Constraints (should be >=0)
    C1 = W_mot - N * Xpmax / pitch  # speed margin
    C2 = T_mot_max - T_mot_scenario  # Torque margin
    
    # --------------------------------------
    # Objective and constraints
    if arg == "objective":
        return objective / 100
    if arg == "objectiveP":
        if ((C1 < 0.0) | (C2 < 0.0)):
            # If constraints are not respected we penalize
            # the objective for contraint less algorithms
            return objective * 1e5
        else:
            return objective / 100
    elif arg == "print":
        print("Objective:")
        print("     Motor mass = %.2f kg" % M_mot)
        print("Design variables:")
        print("     reduction_ratio =  %.2f" % N)
        print("     Rated Motor torque  =  %.2f N.m" % T_mot)
        print("Performances:")
        print("     Max load = %.0f N" % max_load)
        print("Components characteristics:")
        print("     Motor mass = %.2f kg" % M_mot)
        print("     Motor max torque = %.2f N.m" % T_mot_max)
        print("Constraints (should be >= 0):")
        print("     Speed margin: W_mot-reduction_ratio*max_speed/pitch= %.3f" % C1)
        print("     Torque margin:  T_mot_max - T_mot_scenario= %.3f " % C2)
        
    elif arg == "constraints":
        return [C1, C2]
    else:
        raise TypeError(
            "Unknown argument for sizing_code: use 'print', 'objective', 'objectiveP' or 'contraints'"
        )

5.2.4. Optimization with SLSQP algorithm#

We will now use the opmization algorithms of the Scipy package to solve and optimize the configuration. We will first use the SLSQP algorithm without explicit expression of the gradient (Jacobian).

The first step is to give an initial value of optimisation variables:

import numpy as np

# Optimization variables
# Reduction ratio
reduction_ratio_init = 5  # [-]
reduction_ratio_min = 0.1  # [-]
reduction_ratio_max = 10  # [-]

# Oversizing coefficient for multidisciplinary coupling
Tmot_nom = 2  # [-]
Tmot_nom_min = 0.33  # [-] min catalog NK Parvex motor
Tmot_nom_max = 32.8  # [-] max catalog

# Initial values vector for design variables
parameters = np.array((reduction_ratio_init, Tmot_nom))

We can print of the characterisitcs of the problem before optimization with the intitial vector of optimization variables:

# Initial characteristics before optimization
print("-----------------------------------------------")
print("Initial characteristics before optimization :")

sizing_code(parameters, "print")
print("-----------------------------------------------")

-----------------------------------------------
Initial characteristics before optimization :
Objective:
     Motor mass = 2.58 kg
Design variables:
     reduction_ratio =  5.00
     Rated Motor torque  =  2.00 N.m
Performances:
     Max load = 47896 N
Components characteristics:
     Motor mass = 2.58 kg
     Motor max torque = 8.54 N.m
Constraints (should be >= 0):
     Speed margin: W_mot-reduction_ratio*max_speed/pitch= -128.746
     Torque margin:  T_mot_max - T_mot_scenario= -9.077 
-----------------------------------------------

We can see that constraints are not respected.

Then we can solve the problem and print of the optimized solution:

import scipy
import time

# Resolution of the problem

contrainte = lambda x: sizing_code(x, "constraints")
objectif = lambda x: sizing_code(x, "objective")

start = time.time()
result = scipy.optimize.fmin_slsqp(
    func=objectif,
    x0=parameters,
    bounds=[
        (reduction_ratio_min, reduction_ratio_max),
        (Tmot_nom_min, Tmot_nom_max),
    ],
    f_ieqcons=contrainte,
    iter=1000,
    acc=1e-5,
)
end = time.time()

# Final characteristics after optimization
print("-----------------------------------------------")
print("Final characteristics after optimization :")

sizing_code(result, "print")
print("-----------------------------------------------")
print("Calculation time:\n", end - start, "s")

Optimization terminated successfully    (Exit mode 0)
            Current function value: 0.0936881391621254
            Iterations: 17
            Function evaluations: 70
            Gradient evaluations: 17
-----------------------------------------------
Final characteristics after optimization :
Objective:
     Motor mass = 9.37 kg
Design variables:
     reduction_ratio =  2.83
     Rated Motor torque  =  9.00 N.m
Performances:
     Max load = 47896 N
Components characteristics:
     Motor mass = 9.37 kg
     Motor max torque = 38.40 N.m
Constraints (should be >= 0):
     Speed margin: W_mot-reduction_ratio*max_speed/pitch= 0.000
     Torque margin:  T_mot_max - T_mot_scenario= -0.000 
-----------------------------------------------
Calculation time:
 0.1324913501739502 s

/opt/hostedtoolcache/Python/3.9.20/x64/lib/python3.9/site-packages/scipy/optimize/_slsqp_py.py:437: RuntimeWarning: Values in x were outside bounds during a minimize step, clipping to bounds
  fx = wrapped_fun(x)

5.2.5. Optimization with differential evolution algorithm#

We will now use a differential evolution algorithm. Differential Evolution is stochastic in nature (does not use gradient methods) to find the minimium, and can search large areas of candidate space, but often requires larger numbers of function evaluations than conventional gradient based techniques. Differential evolution algorithms don’t manage directly constraints functions. A penalty method replaces the previous objective function.

# Resolution of the problem

objective = lambda x: sizing_code(x, "objectiveP")

start = time.time()
result = scipy.optimize.differential_evolution(
    func=objective,
    bounds=[
        (reduction_ratio_min, reduction_ratio_max),
        (Tmot_nom_min, Tmot_nom_max),
    ],
    tol=1e-5,
)
end = time.time()

# Final characteristics after optimization
print("-----------------------------------------------")
print("Final characteristics after optimization :")

sizing_code(result.x, "print")
print("-----------------------------------------------")
print("Calculation time:\n", end - start, "s")

-----------------------------------------------
Final characteristics after optimization :
Objective:
     Motor mass = 9.37 kg
Design variables:
     reduction_ratio =  2.83
     Rated Motor torque  =  9.00 N.m
Performances:
     Max load = 47896 N
Components characteristics:
     Motor mass = 9.37 kg
     Motor max torque = 38.40 N.m
Constraints (should be >= 0):
     Speed margin: W_mot-reduction_ratio*max_speed/pitch= 0.000
     Torque margin:  T_mot_max - T_mot_scenario= 0.000 
-----------------------------------------------
Calculation time:
 0.13806939125061035 s