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:
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
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