← Back to LocalConstraints1 documentation
LocalConstraints1 - Source Code¶
File: userfiles/SpeedReducer/subsystem1/LocalConstraints1.py
# Copyright (C) The DistributedDesignOptimizer Contributors
# Licensed under the GNU Lesser General Public License v3.0. See LICENSE file for details.
"""Local constraints module for Speed Reducer subsystem 1.
Defines the stress, deflection, and geometric constraints for shaft 1
in the Speed Reducer problem.
"""
import math
from typing import List
from Distributed_Design_Optimizer.subsystem import LocalSubSystemBasis
from Distributed_Design_Optimizer.subsystem.optimization.designproblem import LocalConstraintsInterface
from Distributed_Design_Optimizer.subsystem.tools import ScalerBasis, ScalerConstraint
class LocalConstraints1(LocalConstraintsInterface):
"""Local constraints class for Speed Reducer subsystem 1.
Evaluates shaft 1 stress, deflection, and mounting constraints.
Attributes:
None specific to this class; inherits from LocalConstraintsInterface.
"""
def __init__(self) -> None:
"""Initialize LocalConstraints1 instance."""
pass
def evaluateEqualityLocalConstraints(self, subsystem: LocalSubSystemBasis) -> None:
"""Evaluate equality constraints for subsystem 1.
Args:
subsystem: The local subsystem basis containing state information.
"""
responses: List[float] = subsystem.get_Responses_Unscaled() # unscaled values
scalers: List[ScalerBasis] = subsystem.get_Scalers()
equality_unscaled = []
# append any equality Local constraints to this list using the responses
################################################################
### USER CODE: Equality constraints ###
################################################################
# equality_unscaled.append(responses[...])
# equality_unscaled.append(responses[...])
# scale equality Local constraint evaluation
scl: List[ScalerConstraint] = [scalers[6]]
# equality: List[float] = [scl[i].transform(equality_unscaled[i]) for i in range(len(equality_unscaled))]
# or if no Equality Local constraints exist:
equality = None
# equality Local constraints needs to be a scaled01 quantity
################################################################
### END USER CODE ###
################################################################
subsystem.set_EqualityLocalConstraintsValue(equality)
def evaluate_Jacobian_EqualityLocalConstraints(self, subsystem: LocalSubSystemBasis) -> None:
"""Evaluate the Jacobian of the local equality constraints.
Args:
subsystem: The local subsystem basis containing state information.
"""
return None
def evaluate_Hessians_EqualityLocalConstraints(self, subsystem: LocalSubSystemBasis) -> None:
"""Evaluate the Hessians of the local equality constraints.
Args:
subsystem: The local subsystem basis containing state information.
"""
return None
def evaluateInEqualityLocalConstraints(self, subsystem: LocalSubSystemBasis) -> None:
"""Evaluate inequality constraints for subsystem 1.
Computes 3 shaft-related constraints on stress and geometry.
Args:
subsystem: The local subsystem basis containing state information.
"""
responses: List[float] = subsystem.get_Responses_Unscaled() # unscaled values
scalers: List[ScalerBasis] = subsystem.get_Scalers()
inequality_unscaled = []
# append any inequality Local constraints to this list using the responses
################################################################
### USER CODE: Inequality constraints ###
################################################################
inequality_unscaled.append(((1.93 * (responses[3]**3)) / (responses[1] * responses[2] * (responses[4]**4))) - 1.0)
inequality_unscaled.append((math.sqrt(((745 * responses[3]) / (responses[1] * responses[2]))**2 + 16.9E+6)
/ (110.0 * (responses[4]**3))) - 1.0)
inequality_unscaled.append(((1.5 * responses[4] + 1.9) / responses[3]) - 1.0)
# scale the inequality Local constraint evaluation
scl: List[ScalerConstraint] = scalers[7:10]
inequality: List[float] = [scl[i].transform(inequality_unscaled[i]) for i in range(len(inequality_unscaled))]
# inequality Local constraints needs to be a scaled01 quantity
################################################################
### END USER CODE ###
################################################################
subsystem.set_InequalityLocalConstraintsValue(inequality)
def evaluate_Jacobian_InEqualityLocalConstraints(self, subsystem: LocalSubSystemBasis) -> None:
"""Evaluate the Jacobian of the local inequality constraints.
Args:
subsystem: The local subsystem basis containing state information.
"""
# === Tutorial: scaled <-> unscaled Jacobian (chain rule) ===================
# The optimizer works in SCALED [0, 1] design-variable space, so each Jacobian
# row must be d(scaled constraint) / d(scaled design variables). Every affine
# scaler has a constant slope get_scale() = d(scaled)/d(unscaled). The three
# constraints are written in RESPONSES r = [x1, x2, x3, x4, x6]; the design
# variables are ordered [x4, x6, x1, x2, x3], so we differentiate w.r.t. the
# responses, gather into design-variable order, then convert component-wise:
# dg_s/ds_j = get_scale(constraint) * dg/dx_u[j] / get_scale(dv_scaler_j)
# Return None instead to let the framework fall back to finite differences.
# ===========================================================================
responses: List[float] = subsystem.get_Responses_Unscaled() # unscaled values
scalers: List[ScalerBasis] = subsystem.get_Scalers()
des_var: List[float] = subsystem.get_DesignVariables_Unscaled() # unscaled values
################################################################
### USER CODE: Jacobian of inequality constraints ###
################################################################
# Constraints (unscaled), responses r = [x1, x2, x3, x4, x6]:
# g0 = 1.93*x4^3 / (x2 x3 x6^4) - 1 (monomial)
# g1 = sqrt((745 x4/(x2 x3))^2 + 16.9e6) / (110 x6^3) - 1 (sqrt composite)
# g2 = (1.5 x6 + 1.9)/x4 - 1 (sum of monomials)
dv_to_resp: List[int] = [3, 4, 0, 1, 2]
r: List[float] = list(responses)
n_r: int = len(r)
def mono_grad(coeff: float, p: List[int]) -> List[float]:
M = coeff
for k in range(n_r):
M *= r[k]**p[k]
return [p[k] / r[k] * M for k in range(n_r)]
# g0: monomial in r=[x1,x2,x3,x4,x6]
grad0 = mono_grad(1.93, [0, -1, -1, 3, -4])
# g1 = T * w - 1, with u = 745 x4/(x2 x3), T = sqrt(u^2 + 16.9e6),
# w = (1/110) x6^-3. u = 745 r1^-1 r2^-1 r3^1 (r1=x2, r2=x3, r3=x4).
u = 745.0 * r[3] / (r[1] * r[2])
T = math.sqrt(u**2 + 16.9E+6)
w = (1.0 / 110.0) * r[4]**-3
grad1 = [0.0] * n_r
grad1[1] = w * (u / T) * (-u / r[1]) # dg1/dx2
grad1[2] = w * (u / T) * (-u / r[2]) # dg1/dx3
grad1[3] = w * (u / T) * (u / r[3]) # dg1/dx4
grad1[4] = T * (-3.0 / 110.0) * r[4]**-4 # dg1/dx6
# g2 = 1.5*x6/x4 + 1.9/x4 - 1 (two monomials)
grad2a = mono_grad(1.5, [0, 0, 0, -1, 1])
grad2b = mono_grad(1.9, [0, 0, 0, -1, 0])
grad2 = [grad2a[k] + grad2b[k] for k in range(n_r)]
grads = [grad0, grad1, grad2]
constraint_scalers: List[ScalerBasis] = scalers[7:10]
n_dv: int = len(des_var)
jacobian: List[List[float]] = []
for row in range(3):
scale_c = constraint_scalers[row].get_scale()
jacobian.append([scale_c * grads[row][dv_to_resp[j]] / scalers[j].get_scale()
for j in range(n_dv)])
# jacobian must be a scaled01 quantity (w.r.t. scaled01 design variables)
################################################################
### END USER CODE ###
################################################################
return jacobian
def evaluate_Hessians_InEqualityLocalConstraints(self, subsystem: LocalSubSystemBasis) -> None:
"""Evaluate the Hessians of the local inequality constraints.
Args:
subsystem: The local subsystem basis containing state information.
"""
# === Tutorial: scaled <-> unscaled Hessians (chain rule, 2nd order) ========
# The optimizer works in SCALED [0, 1] space, so each constraint Hessian must
# be d^2(scaled constraint) / d(scaled design vars)^2. Each affine scaler has a
# constant slope get_scale() = d(scaled)/d(unscaled). For the sqrt composite
# g1 = T*w with T = sqrt(u^2 + K), the second derivatives of T follow
# d2T/dx dy = (u_x u_y + u u_xy)/T - u^2 u_x u_y / T^3.
# We build each Hessian in RESPONSE space, gather into design-variable order,
# then convert element-wise:
# d2g_s/ds_a ds_b = get_scale(constraint) * d2g/dx_a dx_b
# / (get_scale(x_a) * get_scale(x_b))
# Return None instead to let the framework fall back to finite differences.
# ===========================================================================
responses: List[float] = subsystem.get_Responses_Unscaled() # unscaled values
scalers: List[ScalerBasis] = subsystem.get_Scalers()
des_var: List[float] = subsystem.get_DesignVariables_Unscaled() # unscaled values
################################################################
### USER CODE: Hessians of inequality constraints ###
################################################################
dv_to_resp: List[int] = [3, 4, 0, 1, 2]
r: List[float] = list(responses)
n_r: int = len(r)
def mono_hess(coeff: float, p: List[int]) -> List[List[float]]:
M = coeff
for k in range(n_r):
M *= r[k]**p[k]
H = [[0.0 for _ in range(n_r)] for _ in range(n_r)]
for a in range(n_r):
for b in range(n_r):
if a == b:
H[a][b] = p[a] * (p[a] - 1) / r[a]**2 * M
else:
H[a][b] = p[a] * p[b] / (r[a] * r[b]) * M
return H
# g0: monomial
H0 = mono_hess(1.93, [0, -1, -1, 3, -4])
# g1 = T * w, T = sqrt(u^2 + K), u = 745 r1^-1 r2^-1 r3^1 (r1=x2,r2=x3,r3=x4),
# w = (1/110) r4^-3 (r4=x6).
K = 16.9E+6
D = 110.0
u = 745.0 * r[3] / (r[1] * r[2])
T = math.sqrt(u**2 + K)
uvars = [1, 2, 3] # response indices appearing in u
pexp = [-1, -1, 1] # corresponding exponents of u
u_x = {vi: pexp[i] / r[vi] * u for i, vi in enumerate(uvars)}
u_xy = {}
for i, vi in enumerate(uvars):
for j, vj in enumerate(uvars):
if vi == vj:
u_xy[(vi, vj)] = pexp[i] * (pexp[i] - 1) / r[vi]**2 * u
else:
u_xy[(vi, vj)] = pexp[i] * pexp[j] / (r[vi] * r[vj]) * u
w = (1.0 / D) * r[4]**-3
wp = (-3.0 / D) * r[4]**-4 # dw/dx6
wpp = (12.0 / D) * r[4]**-5 # d2w/dx6^2
H1 = [[0.0 for _ in range(n_r)] for _ in range(n_r)]
for vi in uvars:
for vj in uvars:
Txy = (u_x[vi] * u_x[vj] + u * u_xy[(vi, vj)]) / T - u**2 * u_x[vi] * u_x[vj] / T**3
H1[vi][vj] = w * Txy
for vi in uvars:
T_i = (u / T) * u_x[vi]
H1[vi][4] = T_i * wp
H1[4][vi] = T_i * wp
H1[4][4] = T * wpp
# g2: sum of two monomials
H2a = mono_hess(1.5, [0, 0, 0, -1, 1])
H2b = mono_hess(1.9, [0, 0, 0, -1, 0])
H2 = [[H2a[a][b] + H2b[a][b] for b in range(n_r)] for a in range(n_r)]
Hs = [H0, H1, H2]
constraint_scalers: List[ScalerBasis] = scalers[7:10]
n_dv: int = len(des_var)
hessians: List[List[List[float]]] = []
for row in range(3):
scale_c = constraint_scalers[row].get_scale()
hessians.append([
[scale_c * Hs[row][dv_to_resp[a]][dv_to_resp[b]] / (scalers[a].get_scale() * scalers[b].get_scale())
for b in range(n_dv)]
for a in range(n_dv)
])
# hessians must be scaled01 quantities (w.r.t. scaled01 design variables)
################################################################
### END USER CODE ###
################################################################
return hessians