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
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:
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
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
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.18/x64/lib/python3.9/site-packages/scipy/optimize/_slsqp_py.py:206, 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)
202 if f_ieqcons:
203 cons += ({'type': 'ineq', 'fun': f_ieqcons, 'jac': fprime_ieqcons,
204 'args': args}, )
--> 206 res = _minimize_slsqp(func, x0, args, jac=fprime, bounds=bounds,
207 constraints=cons, **opts)
208 if full_output:
209 return res['x'], res['fun'], res['nit'], res['status'], res['message']
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/scipy/optimize/_slsqp_py.py:329, in _minimize_slsqp(func, x0, args, jac, bounds, constraints, maxiter, ftol, iprint, disp, eps, callback, finite_diff_rel_step, **unknown_options)
325 # Set the parameters that SLSQP will need
326 # meq, mieq: number of equality and inequality constraints
327 meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
328 for c in cons['eq']]))
--> 329 mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
330 for c in cons['ineq']]))
331 # m = The total number of constraints
332 m = meq + mieq
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/scipy/optimize/_slsqp_py.py:329, in <listcomp>(.0)
325 # Set the parameters that SLSQP will need
326 # meq, mieq: number of equality and inequality constraints
327 meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
328 for c in cons['eq']]))
--> 329 mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
330 for c in cons['ineq']]))
331 # m = The total number of constraints
332 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.18/x64/lib/python3.9/site-packages/scipy/optimize/_differentialevolution.py:392, in differential_evolution(func, bounds, args, strategy, maxiter, popsize, tol, mutation, recombination, seed, callback, disp, polish, init, atol, updating, workers, constraints, x0, integrality, vectorized)
375 # using a context manager means that any created Pool objects are
376 # cleared up.
377 with DifferentialEvolutionSolver(func, bounds, args=args,
378 strategy=strategy,
379 maxiter=maxiter,
(...)
390 integrality=integrality,
391 vectorized=vectorized) as solver:
--> 392 ret = solver.solve()
394 return ret
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/scipy/optimize/_differentialevolution.py:984, in DifferentialEvolutionSolver.solve(self)
979 self.feasible, self.constraint_violation = (
980 self._calculate_population_feasibilities(self.population))
982 # only work out population energies for feasible solutions
983 self.population_energies[self.feasible] = (
--> 984 self._calculate_population_energies(
985 self.population[self.feasible]))
987 self._promote_lowest_energy()
989 # do the optimization.
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/scipy/optimize/_differentialevolution.py:1116, in DifferentialEvolutionSolver._calculate_population_energies(self, population)
1114 parameters_pop = self._scale_parameters(population)
1115 try:
-> 1116 calc_energies = list(
1117 self._mapwrapper(self.func, parameters_pop[0:S])
1118 )
1119 calc_energies = np.squeeze(calc_energies)
1120 except (TypeError, ValueError) as e:
1121 # wrong number of arguments for _mapwrapper
1122 # or wrong length returned from the mapper
File /opt/hostedtoolcache/Python/3.9.18/x64/lib/python3.9/site-packages/scipy/_lib/_util.py:407, in _FunctionWrapper.__call__(self, x)
406 def __call__(self, x):
--> 407 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