Sizing procedure and optimization (Student version)#

Written by Marc Budinger (INSA Toulouse) and Scott Delbecq (ISAE-SUPAERO), Toulouse, France

The objective of this notebook is to learn how to implement a sizing code and use a simple numerical optimization to find the optimal design of the system. The system studied is the TVC EMA of the VEGA launcher.

import numpy as np
import scipy.optimize
from scipy import log10
from math import pi
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[1], line 3
      1 import numpy as np
      2 import scipy.optimize
----> 3 from scipy import log10
      4 from math import pi

ImportError: cannot import name 'log10' from 'scipy' (/opt/hostedtoolcache/Python/3.9.20/x64/lib/python3.9/site-packages/scipy/__init__.py)

Objectives and specifications#

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

The application have to ensure at nozzle level :

  • a max torque \(T_{load}\) of \(48 kNm\) and a max acceleration of \(\dot{\omega}_{max}=811 °/s²\)

  • a max speed \(\omega_{max}\) of 9.24 °/s

  • a max magnitude \(\alpha_{max}\) of 5.7 °

We will give here an example based on a linear actuator with a preselected roller screw (pitch of 10 mm/rev). We assume here, for simplification, the efficiency equal to 70%.

EMA components:

EMA

We first define the specifications and assumptions for the sizing:

# Specifications
angular_magnitude_max = 5.7 * pi / 180  # [rad]
max_dyn_torque = 48e3  # [N.m]
max_speed_rot = 9.24 * pi / 180  # [rad/s]
max_acc_rot = 811 * pi / 180  # [rad/s²]

# Assumptions
pitch = 10e-3 / 2 / pi  # [m/rad]
nu_screw = 0.7  # [-]

# Security coefficient for mechanical components
k_sec = 2
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[2], line 2
      1 # Specifications
----> 2 angular_magnitude_max = 5.7 * pi / 180  # [rad]
      3 max_dyn_torque = 48e3  # [N.m]
      4 max_speed_rot = 9.24 * pi / 180  # [rad/s]

NameError: name 'pi' is not defined

We then define the main characteristics for the components for the scaling laws:

# Motor
T_mot_guess_max_ref = 13.4  # [N.m]
W_mot_max_ref = 754  # [rad/s]
J_mot_ref = 2.9e-4 / 2  # [kg.m²]
M_mot_ref = 3.8  # [kg]

# Rod end
F_rod_max_ref = 183e3  # [N]
M_rod_ref = 1.55  # [kg]
L_rod_ref = 0.061  # [m]

# Screw
M_nut_ref = 2.1  # [kg]
Ml_screw_ref = 9.4  # [kg/m]
D_nut_ref = 0.08  # [m]
L_nut_ref = 0.12 * 0.08 / 0.09  # [m]
F_screw_max_ref = 135e3  # [N]

# Bearing
M_bearing_ref = 5.05  # [kg]
L_bearing_ref = 0.072  # [m]
F_bearing_max_ref = 475e3  # [N]

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.

The design variables of this sizing code are :

  • the reduction ratio of the reducer

  • an oversizing coefficient for the selection of the motor used to tacke an algebraic loop

  • the positions (\(d_1\) and \(d_2\)) of the actuator anchorages

New design variables Mechanical Loads

The objective is the global mass of the actuator.

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

  • the length margin, ie. the global length of the actuator doesn’t exceed the distance between anchorage points

def sizing_code(param, arg="print"):
    #  Design variables
    reduction_ratio = # Reduction ratio of the reducer
    k_oversizing = # Oversizing of the motor (algebraic loop)
    d1 = # Position of the anchorage (cm --> meter)
    d2 = # Position of the anchorage (cm --> meter)

    # --------------------------------------
    # Force, torque and speed calculations

    # Lever arm
    lever_arm = (
        (
            -(-0.9744 * d1 - 1.372)
            * (0.2248 * d1 - 0.3757)
            * ((0.2248 * d1 - 0.3757) ** 2 + (-0.9744 * d1 + d2 - 1.172) ** 2) ** (-0.5)
            + (0.2248 * d1 + 0.9823)
            * ((0.2248 * d1 - 0.3757) ** 2 + (-0.9744 * d1 + d2 - 1.172) ** 2) ** (-0.5)
            * (-0.9744 * d1 + d2 - 1.172)
        )
        ** 2
    ) ** 0.5

    # Length of actuator
    actuator_length = ((0.2248 * d1 - 0.3757) ** 2 + (-0.9744 * d1 + d2 - 1.172) ** 2) ** 0.5

    # Stroke of actuator
    stroke =

    # Load, speed
    max_speed =
    max_load =

    # Torque motor estimation
    T_mot_guess =

    # Max static force (stall force) with 100% efficiency assumption
    max_stall_force =
    max_mech_force =

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

    # Motor
    M_mot =
    J_mot =
    W_mot =

    # Rod end
    M_rod =
    L_rod =

    # Nut
    M_nut =
    M_screw =
    D_nut =
    L_nut =

    # Bearing
    M_bearing =
    L_bearing =

    # --------------------------------------
    # Exact torque calculation with motor inertia
    T_mot_real =

    # --------------------------------------
    # Objectives and constrants calculations

    # Objective = total mass
    objective =

    # Constraints (should be >=0)
    C1 =   # speed margin
    C2 =   # Torque margin
    C3 =   # Length margin

    # --------------------------------------
    # Objective and constraints
    if arg == "objective":
        return objective
    if arg == "objectiveP":
        if () | () | ():
            # 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("     Total mass = %.2f kg" % (M_mot + M_bearing + 2 * M_rod + M_screw + M_nut))
        print("Design variables:")
        print("     reduction_ratio =  %.2f" % reduction_ratio)
        print("     k_oversizing =  %.2f" % k_oversizing)
        print("     d_1 =  %.2f m" % d1)
        print("     d_2 =  %.2f m" % d2)
        print("Performances:")
        print("     Stroke = %.2f m" % stroke)
        print("     Max load = %.0f N" % max_load)
        print("     Stall load = %.0f N" % max_stall_force)
        print("Components characteristics:")
        print("     Lever arm = %.2f m" % lever_arm)
        print("     Actuator length = %.2f m" % actuator_length)
        print("     Motor mass = %.2f kg" % M_mot)
        print("     Max Tem = %.2f N.m" % T_mot_real)
        print("     Rod-end mass = %.2f kg" % (2 * M_rod))
        print("     Rod-end length = %.2f m" % L_rod)
        print("     Screw mass = %.2f kg" % M_screw)
        print("     Nut mass = %.2f kg" % (2 * M_nut))
        print("     Nut length = %.2f m" % L_nut)
        print("     Bearing length = %.2f m" % L_bearing)
        print("Constraints (should be >= 0):")
        print("     Speed margin: W_mot-reduction_ratio*max_speed/pitch= %.3f" % C1)
        print("     Torque margin: T_mot_guess-T_mot_real= %.3f " % C2)
        print("     Length margin: actuator_length-stroke-L_nut-L_bearing-2*L_rod =  %.3f" % C3)
    elif arg == "constraints":
        return [ , , ]
    else:
        raise TypeError(
            "Unknown argument for sizing_code: use 'print', 'objective', 'objectiveP' or 'contraints'"
        )
  Cell In[4], line 3
    reduction_ratio = # Reduction ratio of the reducer
                      ^
SyntaxError: invalid syntax

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:

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

# Oversizing coefficient for multidisciplinary coupling
k_oversizing_init = 1  # [-]
k_oversizing_min = 0.2  # [-]
k_oversizing_max = 5  # [-]

# Anchorage positions
d1_init = 0  # [cm]
d1_min = -80  # [cm]
d1_max = 80  # [cm]

d2_init = 0  # [cm]
d2_min = -20  # [cm]
d2_max = 20  # [cm]

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

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 :
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[6], line 5
      2 print("-----------------------------------------------")
      3 print("Initial characteristics before optimization :")
----> 5 sizing_code(parameters, "print")
      6 print("-----------------------------------------------")

NameError: name 'sizing_code' is not defined

We can see that the consistency constraint that is used to solve the coupling of the motor torque/inertia is not respected.

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

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),
        (k_oversizing_min, k_oversizing_max),
        (d1_min, d1_max),
        (d2_min, d2_max),
    ],
    f_ieqcons=contrainte,
    iter=100,
    acc=1e-8,
)
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")
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 9
      6 objectif = lambda x: sizing_code(x, "objective")
      8 start = time.time()
----> 9 result = scipy.optimize.fmin_slsqp(
     10     func=objectif,
     11     x0=parameters,
     12     bounds=[
     13         (reduction_ratio_min, reduction_ratio_max),
     14         (k_oversizing_min, k_oversizing_max),
     15         (d1_min, d1_max),
     16         (d2_min, d2_max),
     17     ],
     18     f_ieqcons=contrainte,
     19     iter=100,
     20     acc=1e-8,
     21 )
     22 end = time.time()
     24 # Final characteristics after optimization

File /opt/hostedtoolcache/Python/3.9.20/x64/lib/python3.9/site-packages/scipy/optimize/_slsqp_py.py:210, in fmin_slsqp(func, x0, eqcons, f_eqcons, ieqcons, f_ieqcons, bounds, fprime, fprime_eqcons, fprime_ieqcons, args, iter, acc, iprint, disp, full_output, epsilon, callback)
    206 if f_ieqcons:
    207     cons += ({'type': 'ineq', 'fun': f_ieqcons, 'jac': fprime_ieqcons,
    208               'args': args}, )
--> 210 res = _minimize_slsqp(func, x0, args, jac=fprime, bounds=bounds,
    211                       constraints=cons, **opts)
    212 if full_output:
    213     return res['x'], res['fun'], res['nit'], res['status'], res['message']

File /opt/hostedtoolcache/Python/3.9.20/x64/lib/python3.9/site-packages/scipy/optimize/_slsqp_py.py:338, in _minimize_slsqp(func, x0, args, jac, bounds, constraints, maxiter, ftol, iprint, disp, eps, callback, finite_diff_rel_step, **unknown_options)
    334 # Set the parameters that SLSQP will need
    335 # meq, mieq: number of equality and inequality constraints
    336 meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
    337           for c in cons['eq']]))
--> 338 mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
    339            for c in cons['ineq']]))
    340 # m = The total number of constraints
    341 m = meq + mieq

File /opt/hostedtoolcache/Python/3.9.20/x64/lib/python3.9/site-packages/scipy/optimize/_slsqp_py.py:338, in <listcomp>(.0)
    334 # Set the parameters that SLSQP will need
    335 # meq, mieq: number of equality and inequality constraints
    336 meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
    337           for c in cons['eq']]))
--> 338 mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
    339            for c in cons['ineq']]))
    340 # m = The total number of constraints
    341 m = meq + mieq

Cell In[7], line 5, in <lambda>(x)
      1 import time
      3 # Resolution of the problem
----> 5 contrainte = lambda x: sizing_code(x, "constraints")
      6 objectif = lambda x: sizing_code(x, "objective")
      8 start = time.time()

NameError: name 'sizing_code' is not defined

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),
        (k_oversizing_min, k_oversizing_max),
        (d1_min, d1_max),
        (d2_min, d2_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")
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 6
      3 objective = lambda x: sizing_code(x, "objectiveP")
      5 start = time.time()
----> 6 result = scipy.optimize.differential_evolution(
      7     func=objective,
      8     bounds=[
      9         (reduction_ratio_min, reduction_ratio_max),
     10         (k_oversizing_min, k_oversizing_max),
     11         (d1_min, d1_max),
     12         (d2_min, d2_max),
     13     ],
     14     tol=1e-5,
     15 )
     16 end = time.time()
     18 # Final characteristics after optimization

File /opt/hostedtoolcache/Python/3.9.20/x64/lib/python3.9/site-packages/scipy/optimize/_differentialevolution.py:502, in differential_evolution(func, bounds, args, strategy, maxiter, popsize, tol, mutation, recombination, seed, callback, disp, polish, init, atol, updating, workers, constraints, x0, integrality, vectorized)
    485 # using a context manager means that any created Pool objects are
    486 # cleared up.
    487 with DifferentialEvolutionSolver(func, bounds, args=args,
    488                                  strategy=strategy,
    489                                  maxiter=maxiter,
   (...)
    500                                  integrality=integrality,
    501                                  vectorized=vectorized) as solver:
--> 502     ret = solver.solve()
    504 return ret

File /opt/hostedtoolcache/Python/3.9.20/x64/lib/python3.9/site-packages/scipy/optimize/_differentialevolution.py:1155, in DifferentialEvolutionSolver.solve(self)
   1150     self.feasible, self.constraint_violation = (
   1151         self._calculate_population_feasibilities(self.population))
   1153     # only work out population energies for feasible solutions
   1154     self.population_energies[self.feasible] = (
-> 1155         self._calculate_population_energies(
   1156             self.population[self.feasible]))
   1158     self._promote_lowest_energy()
   1160 # do the optimization.

File /opt/hostedtoolcache/Python/3.9.20/x64/lib/python3.9/site-packages/scipy/optimize/_differentialevolution.py:1315, in DifferentialEvolutionSolver._calculate_population_energies(self, population)
   1313 parameters_pop = self._scale_parameters(population)
   1314 try:
-> 1315     calc_energies = list(
   1316         self._mapwrapper(self.func, parameters_pop[0:S])
   1317     )
   1318     calc_energies = np.squeeze(calc_energies)
   1319 except (TypeError, ValueError) as e:
   1320     # wrong number of arguments for _mapwrapper
   1321     # or wrong length returned from the mapper

File /opt/hostedtoolcache/Python/3.9.20/x64/lib/python3.9/site-packages/scipy/_lib/_util.py:441, in _FunctionWrapper.__call__(self, x)
    440 def __call__(self, x):
--> 441     return self.f(x, *self.args)

Cell In[8], line 3, in <lambda>(x)
      1 # Resolution of the problem
----> 3 objective = lambda x: sizing_code(x, "objectiveP")
      5 start = time.time()
      6 result = scipy.optimize.differential_evolution(
      7     func=objective,
      8     bounds=[
   (...)
     14     tol=1e-5,
     15 )

NameError: name 'sizing_code' is not defined