SubSystem Optimization

The Unified Algorithmic Structure defines the following general form for the optimization problem of each LocalSubSystemBasis, together with the Lagrange multiplier associated with each constraint:

\[\begin{aligned} {}^{i}d^{(k,\,l+1)} \leftarrow{} & \argmin_{{}^{i}d\;\in\;{}^{i}\mathcal{D}} \;\; {}^{i}v_f\!\left({}^{i}r\right) + {}^{i}P\!\left({}^{i}d,\; \left\{{}^{i}_{j}u,\;{}^{j}_{i}u\right\}_{j\in{}^{i}N},\; {}^{i}_{C}u,\;{}^{C}_{i}u\right) \\ & \text{s.t.} \;\; {}^{i}v_g\!\left({}^{i}r\right) \leq 0 \quad | \;\; {}^{i}\kappa_g, \\ & \phantom{\text{s.t.}} \;\; {}^{i}v_h\!\left({}^{i}r\right) = 0 \quad | \;\; {}^{i}\kappa_h, \\ & \phantom{\text{s.t.}} \;\; {}^{i}v_{\mathcal{D}}\!\left({}^{i}d\right) \leq 0 \quad | \;\; {}^{i}\kappa_d, \\ & \phantom{\text{s.t.}} \;\; {}^{i}Q^{\leq}\!\left({}^{i}d,\; \left\{{}^{i}_{j}u,\;{}^{j}_{i}u\right\}_{j\in{}^{i}N},\; {}^{i}_{C}u,\;{}^{C}_{i}u\right) \leq 0 \quad | \;\; \kappa_Q^{\leq}, \\ & \phantom{\text{s.t.}} \;\; {}^{i}Q^{=}\!\left({}^{i}d,\; \left\{{}^{i}_{j}u,\;{}^{j}_{i}u\right\}_{j\in{}^{i}N},\; {}^{i}_{C}u,\;{}^{C}_{i}u\right) = 0 \quad | \;\; \kappa_Q^{=}, \end{aligned}\]

where \(\kappa_g\), \(\kappa_h\) and \(\kappa_d\) are the multipliers of the local inequality, local equality and design-variable-bound constraints, and \(\kappa_Q^{\leq}\), \(\kappa_Q^{=}\) are the multipliers of the coordination-specific inequality/equality constraints \(Q^{\leq}\), \(Q^{=}\). All multipliers are gathered at the optimal \({}^{i}\)\(d\) returned by the solver — see Multiplier reconstruction when the solver does not provide them below.

Solving for optimal \({}^{i}\)\(d\)

The Unified Algorithmic Structure Pseudocode-To-Code Traceability indicates that SubSystemBasis.run_IterativeOptimization() solves the above optimization problem. The following procedure details this further.

Highlighted lines are linked to the implementation. Hover or click to see the implementing classes/methods, then click through to the full API documentation.

Procedure SubSystemBasis.run_IterativeOptimization()

optimdata ← OptimizationBasis.callOptimizer(subsystem) ▷ returns optimization results

SolverInterface.execute(subsystem) ▷ dispatched to concrete solver

Retrieve design variable bounds and scaling

 

repeat ▷ solver iterations (objective callback)

Set current design variables

evaluateTotalObjective()

runAnalysis()

AnalysisInterface.evaluateLocalResponses() ▷ user-defined

mapToCouplingParameters()

AnalysisInterface.mapLocalResponsesDesignVariables_to_CouplingParameters() ▷ user-defined

evaluateLocalObjective()

LocalObjectiveInterface.evaluateLocalObjective() ▷ user-defined

evaluateCoordinationObjective() ▷ coordination-specific penalty/augmented terms

Combine: f = LocalObjective + CoordinationObjective

 

evaluateTotalConstraint()

runAnalysis() ▷ same as above

mapToCouplingParameters() ▷ same as above

evaluateLocalConstraints()

LocalConstraintsInterface.evaluateEqualityLocalConstraints() ▷ user-defined

LocalConstraintsInterface.evaluateInEqualityLocalConstraints() ▷ user-defined

evaluateCoordinationEqualityConstraint() ▷ coordination-specific

Combine: ceq = LocalEqConstraints + CoordinationEqConstraints

until solver converged

 

updateSubsystemfromOptimdata(optimdata)

Set optimal design variables, objectives, constraints

mapToCouplingParameters() ▷ final coupling update

runAnalysis() ▷ re-evaluate at optimum

×

The Unified Algorithmic Structure defines the following general form for the optimization problem of each ControllerSubSystemBasis, which is structurally equivalent to the above optimization problem (neglecting \({}^{i}\)\(v_f\), \({}^{i}\)\(v_g\), \({}^{i}\)\(v_h\) and \({}^{i}\)\(D\)):

\[\begin{aligned} \Box \leftarrow{} & \argmin_{\Box} \;\; {}^{C}P\!\left(\Box,\; \left\{{}^{C}_{i}u,\;{}^{i}_{C}u\right\}_{i\in M}\right) \\ & \text{s.t.} \;\; {}^{C}Q^{\leq}\!\left(\Box,\; \left\{{}^{C}_{i}u,\;{}^{i}_{C}u\right\}_{i\in M}\right) \leq 0 \quad | \;\; \kappa_Q^{\leq}, \\ & \phantom{\text{s.t.}} \;\; {}^{C}Q^{=}\!\left(\Box,\; \left\{{}^{C}_{i}u,\;{}^{i}_{C}u\right\}_{i\in M}\right) = 0 \quad | \;\; \kappa_Q^{=}. \end{aligned}\]

Thus the same SubSystemBasis.run_IterativeOptimization() procedure is employed, albeit with coordination-specific definitions of \({}^{C}\)\(d\), \({}^{C}\)\(P\), \({}^{C}\)\(Q^{\leq}\), \({}^{C}\)\(Q^{=}\).

Multiplier reconstruction when the solver does not provide them

The solving workflow described above — SubSystemBasis.run_IterativeOptimization() — is identical for every coordination method and lives entirely in the shared parent classes SubSystemBasis, LocalSubSystemBasis and ControllerSubSystemBasis. Determining the multipliers \(\kappa_g\), \(\kappa_h\), \(\kappa_d\), \(\kappa_Q^{\leq}\), \(\kappa_Q^{=}\) from the KKT conditions, however, is not part of that shared workflow: some coordination methods need these multipliers, while others don't. The KKT determination is therefore only triggered from within the coordination-specific LocalSubSystem<Method>.postprocess_Optimization() (or ControllerSubSystem<Method>.postprocess_Optimization()) — see Where this is used in the workflow below.

Many solvers of the solver/ package do not return these multipliers directly. To make the coordination methods that need them independent of the chosen solver, the subsystem base classes reconstruct them from the first-order optimality (KKT) conditions at the returned primal solution.

The reconstruction logic itself is implemented once in the parent classes SubSystemBasis, LocalSubSystemBasis, and ControllerSubSystemBasis, so that every coordination-specific subsystem that does call it obtains it for free. It relies on the constraint gradients and Jacobians produced by the Derivative Computation, but does not trigger that computation itself: compute_KKT_system_matrix_and_bounds() only reads whatever gradients/Jacobians are already stored in OptimData (e.g. via get_Jacobian_LocalEqualityConstraints()). A coordination method's postprocess_Optimization() must therefore call evaluateAllJacobians() before compute_KKT_multipliers(), as LocalSubSystemALADIN and LocalSubSystemSBDP both do — see Where this is used in the workflow below.

Background: the stationarity condition

At a local optimum \({}^{i}\)\(d\) of the subsystem optimization problem stated above, the KKT stationarity condition states that the gradient of the Lagrangian vanishes:

\[ \nabla_{{}^{i}d} f\left({}^{i}d\right) + \kappa_g^{T} \nabla_{{}^{i}d} {}^{i}v_g\!\left({}^{i}r\right) + \kappa_h^{T} \nabla_{{}^{i}d} {}^{i}v_h\!\left({}^{i}r\right) + \kappa_d^{T} \nabla_{{}^{i}d} {}^{i}v_{\mathcal{D}}\!\left({}^{i}d\right) + \left(\kappa_Q^{\leq}\right)^{T} \nabla_{{}^{i}d} {}^{i}Q^{\leq} + \left(\kappa_Q^{=}\right)^{T} \nabla_{{}^{i}d} {}^{i}Q^{=} = 0, \]

where \(f = {}^{i}\)\(v_f\) \(+\) \({}^{i}\)\(P\) is the subsystem's total objective (local objective plus coordination term, exactly TotalObjectiveValue as assembled by evaluateTotalObjective()), and each gradient term is weighted by the multiplier introduced for that constraint in the optimization problem above. Only active inequality constraints and bounds (and all equality constraints) contribute: an inactive inequality/bound constraint's multiplier is fixed at 0 and drops out of the sum. Writing \(\mathcal{A}\) for this set of active constraints and collecting their gradients row-wise into a Jacobian \(J\) (one row per entry of \(\kappa_g\), \(\kappa_h\), \(\kappa_d\), \(\kappa_Q^{\leq}\), \(\kappa_Q^{=}\) that is active), the stationarity condition becomes the linear system

\[ J^{\top} \lambda = -\,\nabla_{{}^{i}d} f, \]

subject to \(\lambda_a \geq 0\) for the multipliers of active inequality constraints and bounds, and \(\lambda_a\) free (unrestricted sign) for the multipliers of equality constraints. The framework solves this system for \(\lambda\) and then maps its entries back onto \(\kappa_g\), \(\kappa_h\), \(\kappa_d\), \(\kappa_Q^{\leq}\), \(\kappa_Q^{=}\) (or the controller's coordination-only counterparts).

1. Orchestration — compute_KKT_multipliers()

Defined in SubSystemBasis, this is the solver-agnostic driver that assembles and solves the linear system from the previous section, delegating the structure-dependent and fallback work to three helper methods detailed below: it requests the matrix and bounds from compute_KKT_system_matrix_and_bounds() (step 2), forms the right-hand side \(b = -\nabla f\) from the stored gradient of the total objective, solves the system as a linear feasibility problem with a zero cost vector using scipy.optimize.linprog (A_eq \(= J^{\top}\), b_eq \(= b\)), falling back to compute_ApproximateKKT_multipliers() (step 3) if that system is infeasible, and finally hands the solution to decompose_KKT_multipliers() (step 4):

Procedure SubSystemBasis.compute_KKT_multipliers()

result ← compute_KKT_system_matrix_and_bounds() ▷ [constraint_matrix, bounds]

J ← transpose(constraint_matrix)

b ← −gradient_of_total_objective ▷ zeros if no objective gradient

c ← 0 ▷ pure feasibility, no cost

 

result ← linprog(c, A_eq = J, b_eq = b, bounds = bounds)

if status = 0 (success) then all_multipliers ← result.x

else if status = 1 (iteration limit) then raise error

else if status = 2 (infeasible) then warn and

all_multipliers ← compute_ApproximateKKT_multipliers(J, b, bounds)

else raise error ▷ unbounded / unknown status

 

decompose_KKT_multipliers(all_multipliers)

end procedure

An infeasible system (status 2) means either that the primal solution returned by the solver is not exactly optimal, or that the problem does not satisfy constraint qualifications so that no exact multipliers exist. Rather than terminating, the framework prints a warning and falls back to the approximate multipliers of step 3.

2. Assembling the system — compute_KKT_system_matrix_and_bounds()

compute_KKT_system_matrix_and_bounds() is declared abstract in SubSystemInterface and implemented separately in LocalSubSystemBasis and ControllerSubSystemBasis, because local subsystems and the controller carry different constraint sets. It assembles the \(J\) and the per-multiplier bounds (\(\lambda_a \geq 0\) vs. free-sign) of the linear system above: each row of the returned constraint matrix is one active-constraint gradient \(\nabla_{{}^{i}d} c_a\), in the same order as the corresponding entry of \(\kappa_g\), \(\kappa_h\), \(\kappa_d\), \(\kappa_Q^{\leq}\), \(\kappa_Q^{=}\) — a fixed order that both compute_KKT_multipliers() and decompose_KKT_multipliers() rely on:

• LocalSubSystemBasis: ( local equality | coordination equality | local inequality | coordination inequality | lower bounds | upper bounds ), i.e. (\(\kappa_h\)|\(\kappa_Q^{=}\)|\(\kappa_g\)|\(\kappa_Q^{\leq}\)|\(\kappa_d\)(lower) |\(\kappa_d\)(upper)
• ControllerSubSystemBasis: ( coordination equality | coordination inequality | lower bounds | upper bounds ), i.e. (\(\kappa_Q^{=}\)|\(\kappa_Q^{\leq}\)|\(\kappa_d\)(lower) |\(\kappa_d\)(upper) — the controller has no local objective or local constraints.

For each candidate constraint the method decides both the Jacobian row and the multiplier bound:

• Equality constraints contribute their Jacobian row and a free-sign bound (None, None).
• Active inequality constraints and bounds contribute their Jacobian row and a non-negativity bound (0.0, None).
• Inactive inequality constraints and bounds contribute a zero row and the bound (None, None); their multiplier has no impact on the system and is later reported as None.

The activity flags and component Jacobians are read from the subsystem's OptimData (e.g. get_Jacobian_LocalEqualityConstraints(), get_ActiveLowerBounds()), which are populated by the Derivative Computation. The method returns [constraint_matrix, bounds] — constraint_matrix becomes \(J\) (once transposed by step 1 above) and bounds encodes, entry for entry, whether each \(\lambda_a\) is free-sign (equality) or non-negative (active inequality/bound), exactly as required by the stationarity condition above.

3. Fallback — compute_ApproximateKKT_multipliers()

Also defined in SubSystemBasis, this method is called by step 1 only when the exact system is infeasible. It computes the multipliers that minimize the violation of the KKT stationarity condition in the squared \(\ell_2\)-norm,

\[ \min_{\lambda} \; \left\lVert J^{\top}\lambda + \nabla f \right\rVert_2^2 \quad \text{s.t.} \quad \lambda_a \geq 0 \; \text{for active inequality and bound constraints}, \]

by casting it into standard quadratic-program form and solving it with the clarabel solver from the qpsolvers package. The non-negativity requirement of active-constraint multipliers is encoded through the lower bounds; free-sign (equality) multipliers receive a lower bound of \(-\infty\).

4. Distributing the result — decompose_KKT_multipliers()

The solver returns one flat multiplier vector. decompose_KKT_multipliers(), implemented in LocalSubSystemBasis and ControllerSubSystemBasis, walks the vector with an index pointer in exactly the same order used to assemble the matrix and stores each slice on the subsystem's OptimData — for example set_Multipliers_Local_Equality_Constraints() (\(\kappa_h\)), set_Multipliers_Coordination_Equality_Constraints() (\(\kappa_Q^{=}\)), and set_Multipliers_Local_Inequality_Constraints() (\(\kappa_g\)). For inequality and bound constraints, inactive entries are stored as None so that the multiplier list stays aligned with the full constraint list.

Where this is used in the workflow

The computation is triggered lazily and only from within a coordination-specific subsystem class, not from the shared solving workflow: SubSystemBasis.check_MultipliersNotSet() returns True when at least one existing constraint (local bound, local inequality/equality, or coordination equality) is missing its multiplier in OptimData. Only the LocalSubSystem<Method> (or ControllerSubSystem<Method>) classes of coordination methods that actually need the multipliers call the computation from their own postprocess_Optimization(), guarded by this check. Note that compute_KKT_multipliers() itself does not call evaluateAllJacobians() — the caller must run it first so that the Jacobians read by compute_KKT_system_matrix_and_bounds() are up to date:

# LocalSubSystemALADIN.postprocess_Optimization() / LocalSubSystemSBDP.postprocess_Optimization()
self.evaluateAllJacobians()  # populate OptimData's gradients/Jacobians first
self.updateSubsystemfromOptimdata(self.get_OptimData())

if self.check_MultipliersNotSet():
    # Not all multipliers are provided by the solver -> reconstruct them from the KKT conditions.
    self.compute_KKT_multipliers()
self.updateSubsystemfromOptimdata(self.get_OptimData())