# Collaborative Multidisciplinary Design Optimization
with Neural Networks

Jean de Becdelièvre
Stanford University
<EMAIL_ADDRESS>
&Ilan Kroo
Stanford University
<EMAIL_ADDRESS>

###### Abstract

The design of complex engineering systems leads to solving very large
optimization problems involving different disciplines. Strategies allowing
disciplines to optimize in parallel by providing sub-objectives and splitting
the problem into smaller parts, such as Collaborative Optimization, are
promising solutions. However, most of them have slow convergence which reduces
their practical use. Earlier efforts to fasten convergence by learning
surrogate models have not yet succeeded at sufficiently improving the
competitiveness of these strategies. This paper shows that, in the case of
Collaborative Optimization, faster and more reliable convergence can be
obtained by solving an interesting instance of binary classification: on top
of the target label, the training data of one of the two classes contains the
distance to the decision boundary and its derivative. Leveraging this
information, we propose to train a neural network with an asymmetric loss
function, a structure that guarantees Lipshitz continuity, and a
regularization towards respecting basic distance function properties. The
approach is demonstrated on a toy learning example, and then applied to a
multidisciplinary aircraft design problem.

##  1 A Bilevel Architecture for Design Optimization

###  1.1 Introduction

Design optimization of complex engineering systems such as airplanes, robots
or buildings often involves a large number of variables and multiple
disciplines. For instance, an aircraft aerodynamics team must choose the wing
geometry that will efficiently lift the weight of the airplane. This same
weight is impacted by the structures team, who decides of the inner structure
of the wing to sustain the lifting loads. While the disciplines operate mostly
independently, some of their inputs and outputs are critically coupled.

An important body of research on _design architectures_ (Martins and Lambe,
2013) focuses on efficient approaches for solving design problems. Convergence
speed and convenience of implementation are major priorities of the
engineering companies adopting them. Solving design problems as a large
single-level optimization problem has been shown to work well when all
disciplinary calculations can be centralised (Tedford and Martins, 2010).
However, it is often impractical due to the high number of variables and to
the added complexity of combining problems of different nature (continuous and
discrete variables for instance).

Bilevel architectures are specifically conceived to be convenient given the
organization of companies and design groups. An upper system level makes
decisions about the shared variables and assigns optimization sub-objectives
to each disciplinary team at the subspace level. The value of local variables
- such as the inner wing structure - is decided at the subspace level only,
which allows the use of discipline specific optimizers and simulation tools.

Collaborative Optimization (CO)(Braun, 1997; Braun et al., 1996) is the design
architecture that is closest to current systems engineering practices and is
the focus of our work. The system level chooses target values for the shared
variables that minimize the objective function. The sub-objective of each
subspace is to match the targets, or minimize the discrepancy, while
satisfying its disciplinary constraints. Collaborative optimization has been
successfully applied to various large-scale design applications, including
supersonic business jets (Manning, 1999), satellite constellation design
(Budianto and Olds, 2004), internal combustion engines (McAllister and
Simpson, 2003), and bridge structural design under dynamic wing and seismic
constraints (Balling and Rawlings, 2000).

However, the widespread adoption of CO is limited by its slow convergence
speed, which often is an issue with bilevel architectures (Tedford and
Martins, 2010). Mathematical issues when solving the CO system-level problem
with conventional sequential quadratic programming optimizers have been
documented in (Demiguel and Murray, 2000), (Braun et al., 1996) and
(Alexandrov and Robert, 2000), and various adapted versions have been
attempted Roth (2008); Zadeh et al. (2009).

A promising approach is to train surrogate models at predicting whether each
discipline will be able to match the given target value. Using such
representations of the feasible set of each discipline, the system-level
problem can choose better informed target values for the shared variables. The
use of sequentially refined quadratic surrogate models is demonstrated in
Sobieski and Kroo (2000), and several other studies use the mean of a Gaussian
process (GP) regression (Tao et al., 2017). While the earlier work focused on
predicting the discrepancies of each discipline using regression, (Jang et
al., 2005) recognized that the main goal is to classify targets as feasible or
infeasible, and proposes an original mix between a neural network classifier
and a GP regression. Our paper builds upon this work, and shows that the
information obtained from the subspace solution can be used to improve
classification performance.

The rest of the paper is organized as follows. Section 2 shows how CO
reformulates design problems and presents an aircraft design example problem.
Section 3 takes a deeper dive into the structure of the data available to
train a classifier, and proposes to use Lipschitz networks (Anil et al.,
2019). It demonstrates that they generalize better on this task and provides
intuition for it. In Section 4 the example problem of Section 2.2 is solved
with our approach. Finally, Section 5 concludes, discusses the drawbacks of
the current approach and suggests directions for improvements.

##  2 Collaborative Optimization

###  2.1 Mathematical Formulation

Mathematically, the original MDO problem can be written:

$\displaystyle\underset{x_{1:n},z}{\text{minimize}}$ $\displaystyle f(z)$ (1)
subject to $\displaystyle c_{i}(x_{i},z)\leq 0$ $\displaystyle\forall
i\in{1\dots N_{d}}$

where $N_{d}$ is the number of disciplines, $z$ are the shared variables and
$x_{i}$ the local variables of discipline $i$. Typically, the objective $f$ is
easy to compute, but evaluating the constraints $c_{i}$ requires time
consuming disciplinary specific analysis. With CO, we form the _system_ -level
problem that only contains shared variables:

$\displaystyle\underset{z}{\text{ minimize}}$ $\displaystyle f(z)$ (2) subject
to $\displaystyle J_{i}^{*}(z)\leq 0$ $\displaystyle\forall i\in 1\dots N_{d}$

The $(J^{*}_{i})_{i\in{1\dots n}}$ are the values of the _subspace_ -level
optimizations:

$J_{i}^{*}(z)=\min_{{\bar{z}}_{i}\in\\{{\bar{z}}_{i}|\exists
x_{i}\;c_{i}(x_{i},{\bar{z}}_{i})\leq 0\\}}\|{\bar{z}}_{i}-z\|_{2}^{2}$ (3)

${\bar{z}}_{i}$ is a local copy of the global variables. The subspace problem
sets ${\bar{z}}_{i}$ to the feasible point closest to $z$. In practice, the
constraint $J_{i}^{*}(z)\leq 0$ is often replaced by
$J_{i}^{*}(z)\leq\epsilon$, with $\epsilon$ a positive number.

###  2.2 Simple Aircraft Marathon Design Example Problem

We consider the design of an electric powered RC airplane to fly a marathon as
fast as possible. The fuselage, tail, motor and propeller are are given, and
we must design the wing and pick a battery size. There are two disciplines:

  1. 1.

Choosing a wing shape (span $b$ and area $S$) and computing its drag
$D_{wing}$, weight $W_{wing}$ and lift $L$. Making sure $L$ is greater than
the total weight $W=(W_{bat}+W_{fixed}+W_{wing})$. Also, computing the total
aircraft drag $D$.

  2. 2.

Choosing the propulsion system operating conditions (voltage $U$ and RPM) such
that the propeller torque $Q_{p}$ matches the motor torque $Q_{m}$, and
computing the thrust $T$. Making sure $T$ is greater than the total drag
$D=(D_{fixed}+D_{wing})$. Using $V$ and the electric power consumed by the
motor $P_{in}$ to compute $W_{bat}$.

In CO form, the _system_ -level problem only contains the drag $D$ and speed
$V$:

$\displaystyle\min_{D,V,W_{bat}}$ $\displaystyle-V$ subject to $\displaystyle
J^{*}_{aerodynamics}(D,V,W_{bat})=0$ $\displaystyle
J^{*}_{propulsion}(D,V,W_{bat})=0$

where $J^{*}_{aerodynamics}$ and $J^{*}_{propulsion}$ are the optimal values
of the following _subspace_ -problems:
$\displaystyle\min_{\bar{D},\bar{V},\bar{W}_{bat},b,S}$
$\displaystyle\begin{array}[]{c}(\bar{D}-D)^{2}+(\bar{V}-V)^{2}\\\
+(\bar{W}_{bat}-W_{bat})^{2}\end{array}$ subject to $\displaystyle
W_{wing}=\text{Weight}(b,S,L)$ $\displaystyle
L,\bar{D}=\text{Aero}(b,S,\bar{V})$ $\displaystyle L\geq
W_{wing}+W_{fixed}+\bar{W}_{bat}$
$\displaystyle\min_{\bar{D},\bar{V},\bar{W}_{bat},\omega,U}$
$\displaystyle\begin{array}[]{c}(\bar{D}-D)^{2}+(\bar{V}-V)^{2}\\\
+(\bar{W}_{bat}-W_{bat})^{2}\end{array}$ subject to $\displaystyle
Q_{m},P_{in}=\text{Motor}(\omega,U)$ $\displaystyle
Q_{p},T=\text{Propeller}(\omega,\bar{V})$ $\displaystyle
Q_{p}=Q_{m},\;\;T\geq\bar{D}$
$\displaystyle\bar{W}_{bat}=\text{Battery}(\bar{V},P_{in})$

Details about the various models (Weight, Aero, Motor, Battery and Propeller),
as well as a nomenclature can be found in appendix A. This problem is useful
to explain CO, but note that it can easily be solved using single-level
optimization because there is a small number of variables and each discipline
is easily evaluated.

##  3 From the Subspace Problem to Signed Distance Functions

Figure 1: Plots of the square root subspace value function $\sqrt{J^{*}}$ for
an imaginary 1-d problem where the feasible set is $\mathbf{R}^{+}$(left), and
an imaginary 2-d problem where the feasible set is the unit ball (right). The
1-d case also shows the value function $J^{*}$ itself, where the reader can
note the null derivative at the boundary which makes classification ambiguous.
In both cases, $z^{(1)}$ is shown as an example of an infeasible point (along
with its feasible projection ${z^{*}}{(1)}$) while $z^{(2)}$ is shown as a
example of a feasible point.

###  3.1 A Mixed Classification and Regression approach

The subspace optimization problem (Eq. 3), is a projection operation onto the
disciplinary feasible set. Hence, the optimal value of the $i$-th subspace
$J^{*}_{i}(z)$ is the square of the distance between $z$ and closest point to
$z$ inside the feasible set of discipline $i$: it is 0 in the feasible region
and positive outside. Figure 1 shows the aspect of $\sqrt{J^{*}_{i}}$ for a
1-d and 2-d simple examples. The figure also shows hypothetical datapoints
$z_{i}$, as well as the result $z^{*}_{i}$ of the projection.

The main purpose of the $J^{*}_{i}(z)\leq 0$ constraint in the system-level
problem (2) is to indicate the feasibility of $z$. As was noted in Martins and
Lambe (2013) and Tedford and Martins (2010), this constraint is numerically
ambiguous. $J^{*}_{i}(z)\leq\epsilon$ is typically used instead, where the
choice of the threshold $\epsilon$ can negatively impact the feasibility of
the final design. The fact the norm of the gradient of $J^{*}_{i}(z)$ is 0 at
the boundary reinforces the ambiguity. With surrogate models of $J^{*}_{i}$,
the ambiguity becomes even more of a problem. Typical nonlinear regression
approaches do not represent well the large flat areas of the feasible region
and tend to oscillate around zero.

In this work we propose to train a simple feedforward neural network to
classify whether a point $z$ is feasible or not. While exactly regressing on
$J^{*}_{i}$ does not lead to satisfying results, the information about the
distance between each infeasible point and the boundary of the feasible domain
remains very precious. Let $h_{i}$ be a neural network and
$\mathcal{D}=({J_{i}^{*}}^{(j)},z^{(j)})_{j\in{1\dots N}}$ a set of
evaluations of the $i$-th subspace problem (Eq. 3), we propose the following
loss function:

$l(z^{(j)})=\begin{cases}|h_{i}(z^{(j)})-\sqrt{{J_{i}^{*}}^{(j)}}|&\text{if
}z^{(j)}\text{ is infeasible}\\\
\max(h_{i}(z^{(j)}),0)&\text{otherwise}\end{cases}$ (4)

The network is only equal to $\sqrt{J_{i}^{*}}$ in the infeasible region. In
the feasible region, it is trained to be non-positive. The change of sign at
the boundary creates a non-ambiguous classification criterion.

In practice, each subspace solution at $z^{(j)}$ yields a projected point
${z^{*}}^{(j)}$ and the gradient of $J^{*}_{i}(z^{(j)})$. This information can
easily be included in the loss function in equation 4, see appendix B for
details.

Figure 2: Plots of the square root of the subspace value function
$\sqrt{J^{*}}$ as well as the signed distance function of the boundary for an
two imaginary problems. _Top:_ The feasible set is the union of the $[-2,-1]$
and $[1,2]$. _Bottom:_ The feasible set is 2-dimensional and resembles a Pac-
Man character. In this project, we have data about $\sqrt{J^{*}}$, but find it
more efficient for classification to aim at representing the signed distance
function.

###  3.2 Additional Useful Properties of $\sqrt{J_{i}^{*}}$

For any infeasible point $z$, $\sqrt{J^{*}(z)}$ is by definition the distance
between $z$ and the feasible set. Hence, there exist no feasible point in a
ball of radius $\sqrt{J^{*}(z)}$ around $z$. For the trained neural network
$h$ that represents this distance, this property is immediately equivalent to
being 1-Lipshitz:

$\forall(z_{1},z_{2}),\quad\|h(z_{1})-h(z_{2})\|\leq\|z_{1}-z_{2}\|$

Moreover, as can be seen on figure 2, $\sqrt{J^{*}(z)}$ is a continuous
function. But its derivative can be non-continuous if the feasible set is not
convex, which discourages the use of any continuously differentiable
activation function in the neural network $h$.

These properties lead to constraining the type of neural network. In the
remaining examples of this paper, Lipshitz networks designate feedforward
neural networks with orthogonal weight matrices and GroupSort activation
function, as described in Anil et al. (2019). This choice is compatible with
fitting non-continuously differentiable functions and allows to automatically
enforce 1-Lipshitz continuity.

Finally, any function $h$ that represents a distance to a set respects the
following basic properties. For almost every $z$ outside of the set:

$\displaystyle\|\nabla h(z)\|$ $\displaystyle=$ $\displaystyle 1$ (5)
$\displaystyle h\left(z-h(z)\cdot\nabla h(z)\right)$ $\displaystyle=$
$\displaystyle 0$ (6) $\displaystyle\nabla{h\left(x-h(z)\cdot\nabla
h(z)\right)}$ $\displaystyle=$ $\displaystyle\nabla{h(x)}$ (7)

This property is enforced through regularization. At each training iteration,
samples are randomly drawn within the domain, and a term proportional to the
mean squared violation is added to the loss function. Importantly, we apply
this regularization both inside and outside the feasible domain, which
encourages $h$ to represent a _Signed Distance Function_ (SDF) (Osher and
Fedkiw, 2003) of the feasibility boundary, as shown on the right graphs of
figure 2.

###  3.3 Performance comparison as dimension increases

While most fitting approaches do well in low dimension, higher dimensional
problems typically reveal more weaknesses. Building onto the disk example of
figure 1, we let the dimension increase while maintaining the same number of
points. We compare the ability to classify unseen data (as feasible or
unfeasible) of four neural network training approaches:

Name |  Loss function |  Network type
---|---|---
$J^{*}$-fit |  Mean squared error regression on $J$ |  Tanh, feedforward
Classifier |  Hinge loss |  Tanh, feedforward
Hybrid |  Mixed loss (see eqn. 4) |  Tanh, feedforward
SDF |  Mixed loss (see eqn. 4), with
SDF regularization (see eqn. 5) |  Lipshitz network

The results are shown on figure 3, and clearly show that much better
generalization is obtained using the signed distance function approach. The
details of the experiment can be found in appendix C.

Figure 3: Demonstration of the benefit of fitting a signed distance function
to classify disciplinary feasible points given the type of data available in
this work. In 2-d, the top left plots show the train and test data while the
bottom left ones show the contour lines of the learned network in this case.
Fitting $J^{*}$ directly is harder, and also requires choosing an arbitrary
numerical threshold for classification, whereas the SDF only requires a sign
check. On the right, the input dimension is increased while maintaining the
same number of datapoints, demonstrating a much better generalization using
the SDF approach. The experiment is repeated 5 times.

##  4 Solution to the Aircraft Design Problem

The aircraft design problem of section 2.2 is solved using the neural network
classification approach outlined in the previous section. We adopt a very
simple design of experiment:

  1. 1.

Sample points $\left(z^{(j)}\right)_{j\in 1\dots N_{ini}}$ randomly in the
domain and evaluate $J^{*}_{i}(z^{(j)})$ for each discipline $i$.

  2. 2.

Repeat:

    1. (a)

For each discipline $i$, fit a SDF neural network $h_{i}$ following the
methods of section 3.2.

    2. (b)

Find a new candidate point $z$ and add it to the dataset:

$z=\operatorname*{argmin}_{z\in\\{z|\forall i\in 1\dots N_{d}\;h_{i}(z)\leq
0\\}}f(z)$ (8)

If there is no solution to equation 8, we use the point that is closest to
being feasible. While our work focuses on surrogate modeling approaches, more
principled exploration could improve the performance on more complex problems.

Table 1 compares the performance of our approach with a direct application of
CO where the system level problem is solved with sequential quadratic
programming (Braun et al., 1996), as well a with a conventional surrogate
modeling approach using Gaussian process (GP) regression on $J$. Appendix D
gives details about the algorithm and the baselines, and shows visualizations
of the SDF network’s build up.

|  Conventional CO with SQP |  Gaussian Process surrogate model of $J^{i}$ |  Signed Distance approach (ours)
---|---|---|---
No. of system-level iterations |  74 (+/- 28) |  18 (+/- 6) |  6 (+/- 2)
No. of aerodynamic function evaluations |  2161 (+/- 796) |  267 (+/- 63) |  110 (+/- 57)
No. of propulsion function evaluations |  496 (+/- 209) |  94 (+/- 23) |  31 (+/- 8)
Table 1: Solution of the aircraft design problem of section 2.2 using previous
CO approaches and with our approach. The table shows the number of iterations
to 5% of the optimal value. The rounded mean and standard deviation over 20
trials are indicated. The two last column clearly demonstrate the advantage of
using a surrogate model over the original CO architecture. The performance of
our approach comes from the resolution of the ambiguous decision rule
described in section 3.1, as well from a family of function adequate to this
particular problem.

##  5 Conclusion and Future Work

The convergence speed of CO, a distributed multidisciplinary design
optimization architecture is improved by using neural networks to classify
target shared variables as feasible or infeasible. For infeasible points, the
training data contains the distance to the closest feasible point and its
derivative. Leveraging this information by training Lipshitz networks with a
custom asymmetric loss function and proper regularization, we show that more
reliable and faster convergence can be obtained.

An important drawback of the approach presented here is the lack of a
principled exploration strategy. In typical Bayesian optimization (Frazier,
2018; Snoek, 2013; Snoek et al., 2015; Hernández-Lobato et al., 2016), the
exploration-exploitation trade off is carefully taken care of by optimizing an
acquisition function such as expected improvement. This limitation will be
addressed in future work, for instance by building upon the work in (Snoek et
al., 2015).

Future work will also consist in applying our approach to more complex and
larger scale problems. Further extensions could develop multi-fidelity and
multi-objective versions (Huang and Wang, 2009; M. Zadeh and Toropov, 2002),
for which surrogate modeling and Bayesian optimization are already often used
(Picheny, 2015; Rajnarayan et al., 2008; Meliani et al., 2019).

#### Acknowledgment:

This research is funded by the King Abdulaziz City for Science and Technology
through the Center of Excellence for Aeronautics and Astronautics.
https://ceaa.kacst.edu.sa/

## References

  * Martins and Lambe [2013] Joaquim R.R.A. Martins and Andrew B. Lambe.  Multidisciplinary design optimization: A survey of architectures.  _AIAA Journal_ , 51(9):2049–2075, 2013.  ISSN 00011452.  doi: 10.2514/1.J051895.
  * Tedford and Martins [2010] Nathan P. Tedford and Joaquim R.R.A. Martins.  Benchmarking multidisciplinary design optimization algorithms.  _Optimization and Engineering_ , 11(1):159–183, 2010.  ISSN 13894420.  doi: 10.1007/s11081-009-9082-6.
  * Braun [1997] Robert D Braun.  _Collaborative optimization: An architecture for large-scale distributed design._ PhD thesis, Stanford University, 1997.
  * Braun et al. [1996] Robert Braun, Peter Gage, Ilan Kroo, and Ian Sobieski.  Implementation and performance issues in collaborative optimization.  _6th Symposium on Multidisciplinary Analysis and Optimization_ , pages 295–305, 1996.  doi: 10.2514/6.1996-4017.
  * Manning [1999] Valerie Michelle Manning.  _Large-scale design of supersonic aircraft via collaborative optimization_.  PhD thesis, Stanford University, 1999.
  * Budianto and Olds [2004] Irene A Budianto and John R Olds.  Design and deployment of a satellite constellation using collaborative optimization.  _Journal of spacecraft and rockets_ , 41(6):956–963, 2004.
  * McAllister and Simpson [2003] Charles D McAllister and Timothy W Simpson.  Multidisciplinary robust design optimization of an internal combustion engine.  _J. Mech. Des._ , 125(1):124–130, 2003.
  * Balling and Rawlings [2000] R Balling and MR Rawlings.  Collaborative optimization with disciplinary conceptual design.  _Structural and multidisciplinary optimization_ , 20(3):232–241, 2000.
  * Demiguel and Murray [2000] Angel Victor Demiguel and Walter Murray.  An analysis of collaborative optimization methods.  _8th Symposium on Multidisciplinary Analysis and Optimization_ , (c), 2000.  doi: 10.2514/6.2000-4720.
  * Alexandrov and Robert [2000] Natalia M. Alexandrov and Michael Lewis Robert.  Analytical and computational aspects of collaborative optimization.  _NASA Technical Memorandum_ , (210104):1–26, 2000\.  ISSN 04999320.
  * Roth [2008] Brian Douglas Roth.  _Aircraft family design using enhanced collaborative optimization_.  PhD thesis, Stanford University, 2008.
  * Zadeh et al. [2009] Parviz M. Zadeh, Vassili V. Toropov, and Alastair S. Wood.  Metamodel-based collaborative optimization framework.  _Structural and Multidisciplinary Optimization_ , 38(2):103–115, 2009.  ISSN 1615147X.  doi: 10.1007/s00158-008-0286-8.
  * Sobieski and Kroo [2000] I. P. Sobieski and I. M. Kroo.  Collaborative optimization using response surface estimation.  _AIAA journal_ , 38(10):1931–1938, 2000.  ISSN 00011452.  doi: 10.2514/2.847.
  * Tao et al. [2017] Siyu Tao, Kohei Shintani, Ramin Bostanabad, Yu-Chin Chan, Guang Yang, Herb Meingast, and Wei Chen.  Enhanced gaussian process metamodeling and collaborative optimization for vehicle suspension design optimization.  In _International Design Engineering Technical Conferences and Computers and Information in Engineering Conference_ , volume 58134. American Society of Mechanical Engineers, 2017.
  * Jang et al. [2005] B. S. Jang, Y. S. Yang, H. S. Jung, and Y. S. Yeun.  Managing approximation models in collaborative optimization.  _Structural and Multidisciplinary Optimization_ , 30(1):11–26, 2005.  ISSN 1615147X.  doi: 10.1007/s00158-004-0492-y.
  * Anil et al. [2019] Cem Anil, James Lucas, and Roger Grosse.  Sorting out lipschitz function approximation.  In _International Conference on Machine Learning_ , pages 291–301, 2019.
  * Osher and Fedkiw [2003] Stanley Osher and Ronald Fedkiw.  Signed distance functions.  In _Level set methods and dynamic implicit surfaces_ , pages 17–22. Springer, 2003.
  * Frazier [2018] Peter I. Frazier.  A Tutorial on Bayesian Optimization.  2018\.  URL http://arxiv.org/abs/1807.02811.
  * Snoek [2013] Jasper Snoek.  Bayesian Optimization and Semiparametric Models with Applications to Assistive Technology.  2013\.  URL https://tspace.library.utoronto.ca/handle/1807/43732.
  * Snoek et al. [2015] Jasper Snoek, Oren Ripped, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Md Mostofa Ali Patwary, Prabhat, and Ryan P. Adams.  Scalable Bayesian optimization using deep neural networks.  _32nd International Conference on Machine Learning, ICML 2015_ , 3:2161–2170, 2015.
  * Hernández-Lobato et al. [2016] José Miguel Hernández-Lobato, Michael A. Gelbart, Ryan P. Adams, Matthew W. Hoffman, and Zoubin Ghahramani.  A general framework for constrained Bayesian optimization using information-based search.  _Journal of Machine Learning Research_ , 17:1–53, 2016\.  ISSN 15337928.
  * Huang and Wang [2009] Haiyan Huang and Deyu Wang.  A framework of multiobjective collaborative optimization.  In Yong Yuan, Junzhi Cui, and Herbert A. Mang, editors, _Computational Structural Engineering_ , pages 925–933, Dordrecht, 2009. Springer Netherlands.  ISBN 978-90-481-2822-8.
  * M. Zadeh and Toropov [2002] Parviz M. Zadeh and Vassili Toropov.  Multi-fidelity multidisciplinary design optimization based on collaborative optimization framework.  In _9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization_ , page 5504, 2002.
  * Picheny [2015] Victor Picheny.  Multiobjective optimization using gaussian process emulators via stepwise uncertainty reduction.  _Statistics and Computing_ , 25(6):1265–1280, 2015.
  * Rajnarayan et al. [2008] Dev Rajnarayan, Alex Haas, and Ilan Kroo.  A multifidelity gradient-free optimization method and application to aerodynamic design.  _12th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, MAO_ , (September), 2008.  doi: 10.2514/6.2008-6020.
  * Meliani et al. [2019] Mostafa Meliani, Nathalie Bartoli, Thierry Lefebvre, Mohamed-Amine Bouhlel, Joaquim R. R. A. Martins, and Joseph Morlier.  Multi-fidelity efficient global optimization: Methodology and application to airfoil shape design.  2019\.  doi: 10.2514/6.2019-3236.
  * Drela [2007] Mark Drela.  _First-Order DC Electric Motor Model_ , 2007.  downloaded from http://web.mit.edu/drela/Public/web/qprop/motor1_theory.pdf.
  * Brandt et al. [2020] J. B Brandt, R. W. Deters, G. K. Ananda, Dantsker O. D, and M. S. Selig.  _UIUC Propeller Database, Vols 1-3_ , 2020.  downloaded from https://m-selig.ae.illinois.edu/props/propDB.html.
  * Gill et al. [2005] Philip E Gill, Walter Murray, and Michael A Saunders.  Snopt: An sqp algorithm for large-scale constrained optimization.  _SIAM review_ , 47(1):99–131, 2005.
  * Gardner et al. [2018] Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson.  Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration.  In _Advances in Neural Information Processing Systems_ , pages 7576–7586, 2018.

##  Appendix A Details of the Aircraft Design Example Problem

#### Nomenclature

$\bar{X}$ is the subspace-level local copy of the system-level global variable
$X$.

$D_{wing},D_{fixed},D$ | Wing drag, fixed drag (drag of everything but the wing), total drag
---|---
$\eta$ | Propulsion system efficiency
$V$ | Cruise speed
$L$ | Lift
$b,S$ | Wing span and wing area
$W_{wing},W_{bat},W_{fixed},W$ | Wing weight, battery weight, fixed weight (everything but the wing), total weight
$\omega$ | propeller RPM
$T$ | Thrust
$Q_{p},Q_{m}$ | Propeller and motor torque magnitude
$U$ | Motor voltage
$P_{in}$ | Electric power consumed by the motor
$l_{r}$ | Desired Range

#### Propulsion System Model:

The motor voltage $U$ and propeller RPM $\omega$ are the design variables for
this discipline. U is bounded between 0 and 9V, and the $\omega$ between 5000
and 10000.

The motor model is simply a 3 constants model (see for instance [Drela,
2007]). We use a Turnigy D2836/9 950KV Brushless Outrunner (Voltage 7.4V -
14.8V, Max current: 23.2A).

  1. 1.

$K_{v}$ = 950 RPM/V

  2. 2.

$R_{M}$ = 0.07 $\Omega$

  3. 3.

$I_{0}$ = 1.0 A

We use the following model to compute the motor torque and as a function of
the voltage:

$\displaystyle I$ $\displaystyle=(U-\omega/K_{v})/R$ $\displaystyle Q_{m}$
$\displaystyle=\left(I-I_{0}\right)\frac{1}{K}_{q}$

Here we choose $K_{q}$ as equal to $K_{v}$ (converted to units of A/Nm).

The propeller model uses measurements from the UIUC propeller database [Brandt
et al., 2020] for a Graupner 9x5 Slim CAM propeller data. The propeller radius
is $R=0.1143m$. The torque and thrust data for the propeller are fitted using
polynomials:

$\displaystyle Q_{p}$ $\displaystyle=\rho
4/\pi^{3}R^{3}(0.0363\omega^{2}R^{2}+0.0147V\omega R\pi-0.0953V^{2}\pi^{2})$
$\displaystyle C_{T}$ $\displaystyle=0.090-0.0735J-0.1141J^{2}$ $\displaystyle
T$ $\displaystyle=C_{T}*\rho*n^{2}*(2*R)^{4};$ $\displaystyle n$
$\displaystyle=\omega/(2\pi)$ $\displaystyle J$ $\displaystyle=V/n/R$

Finally the battery weight is computed using a battery energy density of
$\nu_{e}=720e3$ J/kg. The power consumed by the propulsive system is
$P_{in}=IU$, and the flighttime for a distance $l_{r}=42000$m is $t=l_{r}/V$,
so we get:

$W_{bat}=\frac{P_{i}nt}{\nu_{e}}=\frac{IUl_{r}}{V\nu_{e}}$

#### Wing Model:

The wing local variables are the span $b$, the wing area $S$, the lift $L$ and
the wing weight $W_{wing}$.

Geometrically the wing is assumed to have a taper ratio $t_{r}=0.75$, a form
factor $k=2.04$, a thickness to chord ratio $t_{c}=0.12$ and a tail area ratio
$t_{t}=1.3$. The area ratio of the airfoil is $k_{airfoil}=0.44$. The fixed
weight of the fuselage, motor and propeller is $W_{f}ixed=37.28N$. The
fuselage has a wetted area $S_{f}=0.18m^{2}$, a body form factor $k_{f}=1.22$
and a length of $l_{f}=0.6m$. We use $g=9.81N/kg$ and $\rho=1.225kg/m^{3}$.

We first use simple aerodynamic theory to compute the lift and the drag of the
wing. We assume $C_{L,max}=1.$ and a span efficiency factor $e=0.8$. The
dynamic viscosity of the air on that day is chosen to be $\nu=1.46e-5m^{2}/s$.
The lift simply is equal to the weight:

$\displaystyle W$ $\displaystyle=W_{fixes}+W_{bat}+W_{wing}$ $\displaystyle L$
$\displaystyle=W$

Then, the induced drag coefficient can be computed:

$\displaystyle C_{L}$ $\displaystyle=2L/(S\rho V^{2})$ $\displaystyle C_{Di}$
$\displaystyle=C_{L}^{2}S/(\pi b^{2}e)$

The parasitic drag is computed assuming fully turbulent flow on the wing and
the fuselage:

$\displaystyle Re_{wing}$ $\displaystyle=\frac{VS}{b\nu}$ $\displaystyle
C_{f,wing}$ $\displaystyle=0.074/Re_{wing}^{0.2}$ $\displaystyle C_{Dp,wing}$
$\displaystyle=(1+2t_{c})C_{f,wing}kt_{t}$ $\displaystyle Re_{fuse}$
$\displaystyle=\frac{Vl_{f}}{\nu}$ $\displaystyle C_{f,fuse}$
$\displaystyle=0.074/Re_{fuse}^{0.2}$ $\displaystyle C_{Dp,fuse}$
$\displaystyle=C_{f,fuse}\frac{S_{f}}{S}k$

Finally we add a penalty $C_{Ds}$ for post-stall flight and complete the final
drag:

$\displaystyle C_{Ds}$ $\displaystyle=0.1\max(0,C_{L}-C_{L,max})^{2}$
$\displaystyle D$ $\displaystyle=\frac{1}{2}\rho
V^{2}S\left(C_{Di}+C_{Dp,wing}+C_{Dp,fuse}+C_{Ds}\right)$

The wing is assumed to be made of a main carbon spar along with styrofoam. The
foam weight is simply related to the wing area. We use
$\rho_{foam}=40kg/m^{3}$

$W_{foam}=g\rho_{foam}\frac{S^{2}}{b}t_{c}k_{airfoil}\\\ $

The gauge of the carbon spar is computed based on stress and deflection, with
a maximum thickness of $\tau_{m}=1.14mm$. The maximum possible stress is
$\sigma_{max}=4.413e9N/m^{2}$, the Young modulus is $E=2.344e11N/m^{2}$ and
the density is $\rho_{carbon}=1380kg/m^{3}$. The radius of the spar is
computed using the wing dimensions: $r_{s}=\frac{S}{4b}t_{c}$. The minimum
carbon thickness to withstand the stress is computed by:

$\displaystyle M_{root}$ $\displaystyle=Lb/8$ $\displaystyle I$
$\displaystyle=\pi r_{s}^{3}\tau_{m}$ $\displaystyle\tau_{stress}$
$\displaystyle=\tau_{m}Lb/8r_{s}/I/\sigma_{max}/0.07$

The minimum carbon thickness to avoid deflections is computed by:

$\displaystyle\delta=Lb^{4}/(64EI)$
$\displaystyle\tau_{defl}=\tau_{m}\frac{2\delta}{b}/0.07$

The actual gauge is the thickest one:
$\tau=\max(\tau_{defl},\tau_{stress},\tau_{m})$. The total mass can then be
computed:

$\displaystyle W_{spar}$ $\displaystyle=2\pi r_{s}\tau b\rho_{carbon}$
$\displaystyle W_{wing}$ $\displaystyle=W_{spar}+W_{foam}$

The optimal value of the problem happens for

D | = 2.15$N$
---|---
V | = 13.71 $m/s$
b | = 2.7574 $m$
S | = 0.397 $m^{2}$
RPM | = 8127
U | = 9$V$

##  Appendix B Including the Subspace Problem Meta-Information in the Loss
Function

As mentioned in the main paper, the subspace optimization problem is a
projection onto the disciplinary feasible set.

$J_{i}^{*}(z)=\min_{{\bar{z}}_{i}\in\\{{\bar{z}}_{i}|\exists
x_{i}\;c_{i}(x_{i},{\bar{z}}_{i})\leq 0\\}}\|{\bar{z}}_{i}-z\|_{2}^{2}$ (9)

[Braun et al., 1996] proves that, for each evaluation $J_{i}^{*}(z)$, we can
compute the gradient of the converged point:

$\nabla J_{i}(z)=(z-{\bar{z}}_{i}^{*})$

where ${\bar{z}}_{i}^{*}$ is the optimal point of the subspace evaluation 9.
We also know from 5 that $J^{i}({\bar{z}}_{i}^{*})=0$, and $\nabla
J^{i}({\bar{z}}_{i}^{*})=\nabla J^{i}(z)$, In cases like section 3 where we
are interested in information about $\sqrt{J^{*}_{i}}$, the gradient
information can be derived easily from this using the chain rule.

A slightly modified version of the loss function of equation 4 is then:

$l(z^{(j)})=\begin{cases}\begin{array}[]{cc}|h_{i}(z^{(j)})-\sqrt{{J_{i}^{*}}^{(j)}}|&+\|\nabla
h_{i}(z^{(j)})-\nabla\sqrt{J_{i}^{*}}^{(j)}\|_{1}+\\\
|h_{i}({({\bar{z}}_{i}^{*})}^{(j)})|&+\|\nabla
h_{i}({({\bar{z}}_{i}^{*})}^{(j)})-\nabla\sqrt{J_{i}^{*}}^{(j)}\|_{1}\end{array}&\text{if
}z^{(j)}\text{ is infeasible}\par\\\
\max(h_{i}(z^{(j)}),0)&\text{otherwise}\end{cases}$ (10)

##  Appendix C Details on the Fits of a Disk of Increasing Dimension

In this experiment, there are 5 different training sets each containing 50
points, while the test set remains the same and contains 500 points. The
experiment is repeated on each train sets in order to obtain error bars on the
reported performance.

All feedforward networks have 3 layers of 8 units each. They are trained with
the Adam Optimizer with $\beta=(0.9,0.999)$ and a $10^{-3}$ learning rate
tuned with gridsearch.

The Lipshitz networks use linear layers orthonormalized after each update
using 15 iterations of Bjorck’s procedure of order 1 with $\beta=0.5$. They
use a full vector sort as activation function (GroupSort with groupsize 1).

##  Appendix D Details on the Solution to the Aircraft Design Problem

###  D.1 Details and Parameter Choices

The system level problem is solved 20 times using each approach, using the
same random seed for each solution approach. The wing and propulsion subspace
problems (see equation 3) are solved using the sequential quadratic
programming solver SNOPT [Gill et al., 2005].

The airspeed $V$ is bound to be between 5 and 15 $m/s$, the battery mass
between $0.1$ and $1$ $kg$, and the drag $D$ between 1 and 6 $N$.

For both surrogate model approaches (Gaussian process and SDF), the initial
number of samples $N_{INI}$ is chosen to be simply 1, which allows to compare
with the direct approach more easily. The choice of a new candidate point is
then done by optimizing the surrogate model, as shown in 8. This optimization
is also performed by SNOPT, with 15 random restarts.

#### Regular CO:

For regular CO, the system level problem is solved directly using SNOPT. A
feasibility threshold of $10^{-4}$ is used to help convergence.

#### Gaussian processes baseline:

At each iteration, a Gaussian Process is fit to the currently available data
for $J_{i}^{*}$, as well as the gradient information. The prior mean is simply
the zero function, and the kernel is a scaled squared-exponential kernel. The
hyper parameters are tuned using type II maximum likelihood minimization and
the Adam Optimizer with $\beta=(0.9,0.999)$ and a $10^{-3}$ learning rate for
200 iterations. The code uses the GpyTorch framework [Gardner et al., 2018].
Again, a feasibility threshold of $10^{-4}$ is required to help convergence.

#### SDF neural networks:

The SDF surrogate models use three layers of 12 units with GroupSort
activation functions (groupsize of 1). The linear layers are orthonormalized
after each update using 15 iterations of Bjorck’s procedure with order 1 and
$\beta=0.5$. They are trained with the Adam optimizer with $\beta=(0.9,0.999)$
and a learning rate that is tuned at each iteration (we use the best network
after out of 10 randomly sampled learning rates between $10^{-5}$ and
$10^{-3}$).