domain: spi phys math info

text:
François Gallard

Charlie Vanaret

Damien Guénot

Vincent Gachelin

Rémi Lafage

Benoit Pauwels

Pierre-Jean Barjhoux

Anne Gazaix

GEMS: A Python Library for Automation of Multidisciplinary Design Optimization Process Generation

This paper describes the GEMS software developed as part of the IRT Saint Exupéry MDA-MDO project 1 for supporting Multidisciplinary Design Optimization (MDO) capabilities. GEMS is a Python library for programing MDO simulation processes, built on top of NumPy 2 , SciPy 3 and Matplotlib 4 .

GEMS aims at pushing forward the limits of automation in simulation process development, with a particular focus on : i) the automatic programming of MDO processes; ii) distributed and multi-level MDO formulations (or MDO architectures); iii) the integration and development of state-of-the-art algorithms for optimization, design of experiments, surrogate models and coupled analyses; iv) the automation of MDO result analysis; v) the deployment of MDO processes in heterogeneous and distributed industrial simulation environments.

Introduction

MDO formulations, define a decomposition strategy to reformulate a design problem into one or multiple optimization problems [START_REF] Martins | Multidisciplinary Design Optimization: A Survey of Architectures[END_REF] . Numerous formulations exist, and their respective performances depend on the addressed problem. For a given problem, their relative performances specially depend on the dimensions of the design variables, coupling variables, and constraints [START_REF] Vanaret | On the Consequences of the "No Free Lunch" Theorem for Optimization on the Choice of an Appropriate MDO Architecture[END_REF] . Besides, industrial applications invoke black box simulation codes that may not be modified or differentiated, or may involve discrete design variables, which restricts the set of potentially useful MDO formulations [START_REF] Gazaix | Towards the Industrialization of New MDO Methodologies and Tools for Aircraft Design[END_REF] . In addition, the deployment of MDO at industrial scale implies implementing numerous MDO processes, which can be very expensive if not supported by appropriate software.

The GEMS library was developed to tackle these industrial constraints. MDO processes should be generated automatically with a minimal amount of programming effort, based on a catalog of available MDO formulations which may be adapted to dedicated use cases.

In the aeronautical field, disciplinary optimization software, such as aerodynamic shape optimization or structural sizing, have typically more advanced capabilities than MDO capabilities involving both highfidelity aerodynamics and structure simulations in a single process. In the past decades, considerable efforts were committed to the industrialization of such disciplinary problems solvers [START_REF] Johnson | Thirty years of development and application of CFD at Boeing Commercial Airplanes, Seattle[END_REF][START_REF] Meaux | Viscous aerodynamic shape optimization based on the discrete adjoint state for 3d industrial configurations[END_REF][START_REF] Grihon | Numerical Optimization applied to structure sizing at AIRBUS: A multi-step process[END_REF] . The methodologies used in these solvers take advantage of the disciplinary optimization problem structure, such as the linearity of the mass objective and the large amount of non-linear constraints in structural sizing problems, while the aerodynamic optimization problems have a non-linear objective function, the drag, with a small number of non-linear constraints such as the lift or only bounds on design variables. This makes them extremely efficient, but specialized. If one wants to reuse these capabilities in a MDO process, such a process shall involve sub-optimization problems and therefore the selected MDO formulation shall be distributed. Multilevel MDO formulations decompose the MDO problem into multiple optimization problems. For instance in the BLISS98 [START_REF] Sobieszczanski-Sobieski | Bi-Level Integrated System Synthesis (BLISS)[END_REF] , the design variables that are shared by multiple disciplines are the unknowns of a system-level optimization, while the design variables that are only directly impacting a single discipline are the unknown of so-called disciplinary optimizations. These can therefore make direct use of existing disciplinary optimization software as opposed to monolithic MDO formulations, such as MDF, which involve a single optimization problem with coupled disciplines. Another advantage of this approach is that the MDO formulations are independent of the disciplinary tools, so a single implementation can be validated and shared for multiple use cases.

For these reasons, GEMS was developed with multi-level MDO formulations as an important use case, with the possibility to reuse existing black-box optimization software. This is a substantial difference with NASA's OpenMDAO [START_REF] Gray | OpenMDAO: An open source framework for multidisciplinary analysis and optimization[END_REF] , which mainly targets MDF and IDF since the rewriting of the code for version 1.0.0 [START_REF] Gray | Automatic evaluation of multidisciplinary derivatives using a graph-based problem formulation in OpenMDAO[END_REF] . The focus is made in OpenMDAO on getting the maximal performance out of these formulations, taking advantage of parallel algorithms and the graph structure of the process. The KADMOS software [START_REF] Van Gent | Composing MDAO symphonies: graph-based generation and manipulation of large multidisciplinary systems[END_REF] , has MDO formulations capabilities, and is also mainly focused on MDF and IDF formulations as of 2017 in version 0.7. In comparison, GEMS is a generic library that supports both monolithic and multi-level MDO formulations. Another difference with the KADMOS software is that GEMS does not delegate the execution of the MDO process to a workflow engine but achieves itself the workflow and dataflow.

In Section 2, we propose a review of the existing paradigms for programming MDO processes, and conclude on the need for a new approach, implemented in GEMS. In Section 3, we describe GEMS' software architecture. In Section 4, we list some of the available features, among which the MDO formulations, the MDO results visualisations, and the MDA (Multidisciplinary Design Analysis) algorithms with coupled derivatives.

MDO processes programing : state of the art and a new paradigm

The basics of simulation process programming

In order to build an automated simulation process, the elementary bricks of the simulation processes shall be interfaced with a workflow manager, itself responsible for the overall execution and exchange of data. These interfaces are achieved through wrappers, which are a piece of software dedicated to the execution of other programs. Each program is encapsulated using a dedicated interface. A wrapper is an abstraction that is used to define input data, output data and an execution of the interfaced software independently of its own implementation (file formats and data structures for instance). The concept is used in all process integration software (OpenMDAO from NASA [START_REF] Gray | OpenMDAO: An open source framework for multidisciplinary analysis and optimization[END_REF] , ModelCenter TM from Phoenix Integration 14 or Isight TM from Dassault Systèmes [START_REF] Systèmes | SIMULIA Isight 5.7 User's Guide[END_REF] ). It is a prerequisite to build any process.

The motivations of programming simulation processes using wrappers and a worfklow engine, compared to scripting it directly with a programming language, are the following.

• Easily create and maintain many similar processes using the same software wrappers.

• Make the process less dependent on the interfaced software version changes: input and output modification, internal methods implementation.

• Develop automated checks of the consistency of the process to anticipate errors before execution.

• Develop common preparation and analysis tools for multiple processes.

Data-driven versus workflow-driven

Several paradigms can be used to program a simulation process from wrapped software. The simulation process is defined by the execution sequence of the wrapped software, called "the workflow", and the data exchanges between them, called "the dataflow". There are many possible classifications of workflow engines: we propose here one particular view which is well suited to present the GEMS concepts in the next sections. Currently existing workflow engines can be classified into three paradigms for programming the simulation processes.

1. Full process description: the process developer programs both the data exchanges between the wrapped software, and the execution sequence. The process is neither driven by the data nor the workflow, since both are inputs of the workflow engine. This is the case of ModelCenter TM from Phoenix Integration 14 , where the process developer shall first define the order of execution, and then the data exchanges between the boxes.

2. Data-driven: the process builder describes only the data exchanges of the process, and the workflow engine automatically generates the workflow. The process is then "data-driven" because the data exchanges drives the whole process. This is typically the case of MATLAB Simulink TM processes.

3. Workflow-driven: the process builder describes only the sequence of execution of the process, and the workflow engine automatically generates the dataflow. The process is then "workflow-driven" because the execution sequence drives the whole process. This is typically the case of the WORMS workflow engine, developed at Airbus by Meaux et al. for aerodynamic shape optimization [START_REF] Meaux | Viscous aerodynamic shape optimization based on the discrete adjoint state for 3d industrial configurations[END_REF][START_REF] Meaux | Amélioration d'une chaine d'optimisation numérique de formes aérodynamiques en vue du traitement de configurations complexes en écoulements visqueux[END_REF][START_REF] Gallard | Optimisation de forme d'un avion pour sa performance sur une mission[END_REF] , and later eWORMS, based on Eclipse RCP.

In the following paragraphs, we examine the pro and cons of data-driven versus workflow-driven engines in the perspective of MDO for aeronautical test cases.

Full process description

Figure 1 shows how the execution sequence and the data exchanges of the process are defined graphically by the builder of a ModelCenter process. The two wrapped software A and B are executed in the order A then B, as shown in Figure 1a. Then, the output "o2" of A, is connected to the input "o1" of B, as shown in Figure 1b. In this case, the process builder programs both the workflow and the dataflow of the process. The workflow engine is therefore not automating a part of the process programming, but the process execution only. However, the builder has full direct control over the workflow and dataflow, which can be useful when the interfaced software do not use consistent names for their inputs and outputs, nor follow a standardized or generic process.

The advantages of workflow engines based on the full process description are:

• the process builder can fully customize the data exchanges and execution sequence to specific cases;

• there are no rules to be followed on the data names, nor generic process to be defined.

The main weaknesses of the approach are the following:

• When the interfaced software have many inputs and outputs, the process becomes very hard to build.

• The workflow is very sensitive to inputs and outputs changes. If, during a version change of a component, a new input appears, the whole process may be impacted. The process has a weak tolerance to inputs and outputs changes.

• The workflow is impacted by a change of interfaced software, changing one of them impacts the workflow execution sequence.

Data-driven

The data-driven paradigm is inspired from electrical wiring and discrete time simulation to automate the execution sequence of the process: the wrapped software are executed every time their inputs are ready. Figure 2 illustrates the data exchanges specified by a process builder for a typical simulation process built with a data-driven workflow engine such as CFD or CSM. The steps of the process are:

1. building the geometry, 2. generating the mesh from the geometry, 3. solving the equations of interests (RANS equations for instance in CFD) on the generated mesh, 4. post-processing the simulation (computing lift and drag for instance).

The data-driven way to describe this process is to describe all the exchanges between the wrapped software in terms of data. To illustrate the process execution in this particular case, when the user supplies a new set of parameters to the geometry, the geometry program generates a file that is passed to the meshing program, which triggers its execution, and produces a mesh. This mesh is then passed to the simulation program, along with geometrical information, which triggers the simulation. Finally, the solution, the mesh and the geometrical information are passed to the post-processing tool to compute quantities of interest. In data-driven workflow engines, the data drive the process: the builder specifies the data exchanges and the workflow is automated by the engine. Therefore, the workflow engine has specific rules to decide when the next wrapped program should be executed, that is when its inputs are ready. Each workflow engine has its own rules depending on the application. For signal processing for instance, fixed time steps may be used and all executions may be triggered simultaneously for each step.

The advantages of data-driven workflow engines are:

• the workflow is automated, the process developer does not have to care about the overall execution sequence;

• focus is on the data;

• it is a well suited approach for signal processing and electrical system designs, which have many components and low-dimensional data.

The main weaknesses of the approach are the following:

• When the interfaced software have many inputs and outputs, the process becomes very hard to build.

• The workflow has a weak tolerance to input and output changes. If a new input appears during a version change of a component, the whole process may be impacted.

Workflow-driven

In high-fidelity simulations, the inputs and outputs of the program are numerous, a CFD code for instance has typically hundreds of inputs and options. To overcome the limitations of data-driven paradigm and achieve a partial automation of the process development, workflow-driven engines have been developed.

The main idea of workflow-driven engines is that, in some processes such as CFD simulations, the execution order is constant. The geometry will always be generated before the mesh, which will then be used by the CFD solver, and the postprocessing will always come last. The predefined execution order, along with rules on the data exchanges, defines the dataflow. The workflow engine is in charge of managing the data exchanges accordingly. Namespaces as well as data privacy management between programs are important points to build complex processes.

Figure 3 illustrates the execution sequence specified by a process builder for a typical simulation process built with a workflow-driven engine. The advantages of workflow-driven workflow engines are:

• the dataflow is automated, the process developer does not have to care about the overall data exchanges, which is good for processes with a lot of data;

• focus is on the process execution order, which is well suited for some classes of processes such as CFD;

• it is a well suited paradigm for processes with a constant execution order, or programs whose inputs and outputs change over versions.

The main weaknesses of the approach are the following:

• The user has difficulties to have a clear representation of the process. Most people are used to boxes and wires, that is a data-driven way to display processes, which is not suited for workflow-driven engines, where the process is a tree of execution sequences.

• The process is sensitive to change in the interfaced software. For instance, if one wants to replace a simulation post processing tool by another, all processes that use the former will have to be reprogrammed.

Dataflow and workflow automation needs for MDO

From the previous sections, we can see that there is a duality between the dataflow and the workflow in the current way of programming simulation processes. In either case, the dataflow or/and the workflow have to be programmed. In high-fidelity MDO for aerospace applications we have disciplines with many inputs and outputs, which change over time and use cases. Therefore the process programmer can hardly care about all the data. In addition, for given disciplines, many different problems and associated processes are built in an industrial-scale development. Besides, MDO formulations exhibit different performances when the use case changes, since there is no "free lunch" for MDO formulations [START_REF] Vanaret | On the Consequences of the "No Free Lunch" Theorem for Optimization on the Choice of an Appropriate MDO Architecture[END_REF] . For given disciplines, changing the formulation completely changes the dataflow and the workflow. Figure 4 shows some use cases addressed in the MDA-MDO project 1 . For each use case, multiple methods may be used, such as DOEs, RSMs, or MDO formulations, and for each methods, different algorithms may also be selected. Besides, the automation of the process shall be implemented in a workflow engine, and we display here three possible options (Isight, ModelCenter, WORMS). There are 99 paths in this graph: each one represents a possible process, i.e. a combination of a use case, a formulation, an algorithm and a workflow engine. In addition, each of these items is implemented in a separate software package, which evolves in time. For three different versions of each software, this leads to 3 99 ≈ 10 47 potential processes. This is a major issue for maintenance and industrial deployment of MDO methods, which illustrates why MDO process programming must be automated and reconfigurable. As a conclusion, we need automated dataflows as in workflow-driven engines and automated workflow management as in data-driven engines. Consequently, new paradigms to implement MDO simulation processes are required to overcome the current limitations.

GEMS' solution to address both dataflow and workflow automation

Some of the limitations of workflow-driven engines are solved by using the object oriented programming paradigm. If wrappers are classes, one can implement subclasses by inheritance, which then allows to substitute a subclass by another using polymorphism. This allows for instance to change one of the simulation tool by another with identical purpose in a process, but not to reconfigure the whole process with a completely different execution sequence. One cannot reconfigure in this way a MDF process, based on a single optimization calling a MDA, into a bi-level MDO process, based on multiple disciplinary optimizations.

With the objective of developing reconfigurable MDO processes in mind, generic processes have been developed in GEMS, which can be used for any wrapped program. MDAs (Multi-Disciplinary Analysis) for instance are generic processes, so are software chains that represent a linear execution sequence of a list of disciplines. Once instanciated with a list of disciplines as attributes, these generic process can be executed to achieve a workflow. They are programmed in a workflow-driven way to take advantage of this paradigm.

A second limitation in the current way of building the processes is the link with optimization algorithms. The definitions of the objective function and constraints depend on each test case, since the quantities of interest are actually computed by different programs. Either in workflow or data-driven engines, the optimization loops are programmed within the rest of the process. Therefore, the formulation is problemspecific, which is an issue for reconfiguring the process to use a different MDO formulation. Besides, current workflow engines have succeeded to support disciplinary optimization processes, but they are unadapted to multi-level MDO formulations, which involve multiple nested optimization loops. And above all, they are unable to support a complete process reconfiguration such as going from a monolithic formulation to a multi-level formulation.

MDO formulations are generic strategies that generate optimization problems on top of generic processes. Objective functions and constraints can be computed from the processes outputs. This is a generic approach that can be implemented independently of the test case. The dataflow is handled in a part by the optimization algorithm and in a part by the generic processes, as a consequence of the selected MDO formulation. The same applies for multi-level formulations such as BLISS98 [START_REF] Sobieszczanski-Sobieski | Bi-Level Integrated System Synthesis (BLISS)[END_REF] , since the decomposition strategy of the MDO problem into a system level and disciplinary optimization problems is also generic.

High-fidelity MDO processes take weeks to run. The traditional trial-and-error way to debug such a process is unrealistic. The GEMS way of building processes allows a multi-step for the validation and debugging of MDO processes:

• the generic processes and MDO formulations are independent of the test case, so they can be validated on well known test cases;

• the wrappers are independent of the processes, so the integration of the simulation programs can be validated separately;

• the performance of algorithms and formulations can be tuned on substitution problems thanks to the scalable methodology by Vanaret et al. [START_REF] Vanaret | On the Consequences of the "No Free Lunch" Theorem for Optimization on the Choice of an Appropriate MDO Architecture[END_REF] 3. Software architecture

General principles

In this section, some parts of GEMS' design are discussed, with focus on some key classes. GEMS has more than 100 classes, therefore the full class diagram is not shown here. The package has an intermediate size of 26,000 lines of Python code. Python was chosen because the main purpose of GEMS is to integrate multiple software, and Python is known to be efficient to glue together existing components in heterogeneous languages when compared to C++ and Java [START_REF]Glue It All Together With Python[END_REF] . Besides, GEMS also contains numerical algorithms manipulating matrices and this can be programmed easily using he NumPy 2 , SciPy 3 libraries. Many of the basic numerical algorithms for scientific computing (linear and non-linear systems, optimisation, signal processing) are provided in these libraries. Finally, optimization data plots can be easily programmed from the NumPy arrays produced by the scenarios execution using the Matplotlib library 4 .

Main classes

The Scenario (DOE or MDO) is composed of disciplines, and the MDO formulation. The MDO formulation creates an optimization problem from the list of disciplines and a design space. This may require the intermediate creation of generic sub-processes such as MDAs, or other MDO scenarios. All generic processes inherit from disciplines and are composed of a list of disciplines. An optimization problem is defined from functions generated by the FunctionGenerator, on top of the generic processes or directly the disciplines as in the case of the IDF formulation. The drivers (optimization algorithm or DOE) solve or sample the optimization problem, and when started, trigger the execution of the process.

The high-level classes that are key in the architecture are the following.

• MDOScenario builds the process from a set of inputs, several disciplines and a formulation. It is one of the main interface classes for the MDO user. The MDOScenario triggers the overall optimization process when its execute() method is called. Through this class, the user provides an initial solution, user constraints and the design space. The user may also generate visualizations of the scenario execution, such as convergence plots of the algorithm.

• DOEScenario is similar to the MDOScenario, but is dedicated to DOE-based studies. Both MDOScenario and DOEScenario inherit from the Scenario class that defines common features.

• MDODiscipline represents a wrapped simulation program or a chain of wrapped software. It can either be a link to a discipline interfaced with a workflow engine, or can be inherited to interface a simulation program directly. Its inputs and outputs are represented in a Grammar.

• MDOFormulation describes the MDO formulation (e.g. MDF and IDF) used to solve the MDO problem.

• OptimizationProblem describes the mathematical functions of the optimization problem (objective function and constraints), along with the design variables. It is generated by the MDOFormulation, and solved by the optimization algorithm. It has an internal database that stores the calls to its functions by the optimization algorithm to avoid duplicate computations. It can be stored on disk and analyzed by the post-processing methods.

• DesignSpace is an attribute of the OptimizationProblem that describes the design variables, their bounds, their type (floating point or integer), and current value. This object can be read from a file.

Two low-level classes are crucial to understand the basic principles of GEMS:

• MDOFunction instances (objective function and possible constraints) are generated by the MDOFormulation.

Depending on the formulation, constraints may be generated as well (e.g. consistency constraints in IDF).

• MDOFunctionGenerator handles the creation of MDOFunctions for a given MDODiscipline. It is a key class for the development of MDO formulations because these functions are passed to the optimizers, which drive the execution, and their values and gradients are computed from the disciplines. Therefore, the functions are a bridge between the MDO formulations and the generic processes part of GEMS. Figure 5 shows the composite design pattern used for the MDOScenario, which is a subclass of MDODiscipline and is also composed of a list of MDODisciplines. It is a key point for multilevel MDO formulations because sub-MDOScenarios can be embedded into a system-level MDOScenario as its MDODisciplines.

Figure 7 shows the abstraction used to define DOEScenario and MDOScenario that both inherit from the main class Scenario. A DOEScenario can also be substituted with an MDOScenario, so trade-off processes can be built using a bi-level MDO formulation that mixes design of experiments and disciplinary suboptimizations.

Figure 8 shows a part of the classes hierarchy for the generic processes. All generic processes are subclasses of the MDODiscipline class. This composite design pattern is the key for building complex processes by blocks. The Gauss-Seidel algorithm for instance can be used to solve an MDA, and is built using a sequential execution ("chained") of the disciplines by the MDOChain class. The MDA is also a subclass of MDODiscipline, so can be embedded in any other generic process as a substitute of a MDODiscipline. Advanced MDAs as shown in Section 4.7.1 are composed of a MDOChain, taking parallel execution, MDAs and single MDODiscipline as its own input disciplines, in order to achieve the graph of execution plotted in Figure 22.

This building blocks design is also exploited for the automated derivatives computation. The computation of coupled derivatives is implemented once in the MDA base class and then shared by all MDA sub-classes. Besides, the classical chain rule for composing derivatives is available in the MDOChain class. Any process that is built on top of these classes will automatically benefit from these capabilities. Advanced MDAs, implemented MDAChain, can therefore either compute the derivatives of the overall coupled problem through its base MDA class, or by chain ruling the coupled derivatives of the sub-MDAs through the MDOChain.

Execution sequence

Figure 9 shows the execution sequence diagram of an MDF scenario. The scenario triggers the overall execution through the resolution of an optimization problem. The generated objective and constraints are passed to the optimizer, which calls them iteratively. At each call, a MDA is executed. A caching mechanism at the optimization problem level is also implemented to avoid duplicate executions, and is not illustrated here. Then, the algorithms may also call the gradient of the objective function and constraints, which triggers the "linearize" method of the MDA, which also needs the Jacobian matrices of the disciplines in order to compute the coupled derivatives. These Jacobian matrices are computed by the "linearize" methods of the disciplines.

A similar scheme applies to multi-level scenarios. At system level, the MDA of Figure 9 is replaced by multiple MDO scenarios. Figure 10 shows how GEMS (the MDO formulations engine) is typically interfaced in an MDO platform. The platform provides standardized interfaces for disciplinary workflow engines that run on HPCs. This allows to build a multi-machine distributed MDO capability. GEMS has standard interfaces to optimization and DOE algorithms. A proof of concept of a distributed MDO process running on multiple machines under Windows or Linux and with different workflow engines has been performed 1 . The target platform for such a component-based architecture in the project is Eclipse RCP, [START_REF]Eclipse Rich Client Platform[END_REF] which allows graphical integration of multiple plugins. Such a graphical user interface is being developed on top of GEMS to provide user inputs, configure and monitor a MDO scenario.

Features overview

MDO formulations

The following MDO formulations have been developed with GEMS:

• Multidisciplinary Design Feasible (MDF).

• Individual Design Feasible (IDF).

• Disciplinary optimization for weakly coupled problems.

• Two variants of bi-level formulation using disciplinary design scenario 1 , derived from ONERA's research on high-fidelity MDO for aeronautical applications 20 .

• A bi-level formulation based on distributed MDF scenarios 1 .

• A bi-level formulation based on asymmetric distributed MDF scenarios 1 .

• A distributed MDO formulation with overlap on the shared design variables, based on a domain decomposition method with synchronization 21 .

• A distributed formulation dedicated to structural optimization with sizing and categorical variables based on the branch and bound algorithm [START_REF] Barjhoux | Mixed Variable Structural Optimization: toward an Efficient Hybrid Algorithm[END_REF] .

Once the MDO process is built (however not executed), GEMS can generate the XDSM diagram 23 associated to the selected MDO formulation using ONERA's XDSMjs library 24 . This feature is useful to the user to represent himself a simplified view of the process that is generated by GEMS. It also allows to check that all the connections (the dataflow and the workflow) are achieved properly before execution, and to debug the process.

GEMS scenario execution

Appendix A and Appendix B show the output log of two MDO scenarios executed with GEMS on Sobieski's SSBJ 10 design optimization problem. The initial design space is provided, as well as the generated optimization problem. Additional consistency constraints are generated by the IDF formulation. Note that the range objective function (y 4 ) and the design constraints (g 1 , g 2 , g 3 ) have different dependencies due to the different formulations. During the execution, a progress bar indicates the optimization progress as well as the number of iterations per second. The best feasible point is given as the solution of the problem. Finally the execution metrics can be printed, such as the total spent time and number of calls to each discipline.

Appendix C shows the first system iteration of the first bi-level process described by Gazaix et al. [START_REF] Gazaix | Towards the Industrialization of New MDO Methodologies and Tools for Aircraft Design[END_REF] whose XDSM [START_REF] Lambe | Extensions to the design structure matrix for the description of multidisciplinary design, analysis, and optimization processes[END_REF] is also displayed in Figure 11. One can note the system optimization problem that handles the shared design variables, and the disciplinary optimizations. MDA histories are not shown here.

Interfaces to a range of DOE and optimization algorithms

Once GEMS has generated an optimization problem, its resolution is delegated to an optimization algorithm. A set of optimization and DOE algorithms are therefore interfaced with GEMS. The high performance of an optimization algorithm measured on a set of optimization problems necessarily comes with a lower relative performance on other sets of optimization problems, which is known as the "No free lunch theorem for optimization" [START_REF] Wolpert | No free lunch theorems for optimization[END_REF] . This is particularly important in MDO since changing the coupling or design variables dimensions generates different optimization problems, which has an impact on the performance of the algorithms and in the end of the MDO process [START_REF] Vanaret | On the Consequences of the "No Free Lunch" Theorem for Optimization on the Choice of an Appropriate MDO Architecture[END_REF] . It is therefore key to allow the user to interface easily a dedicated optimization algorithm for solving its particular problem. GEMS provides a standardized interface to this aim. Currently 16 optimization algorithms are available from the NLopt 26 , OpenOpt 27 , SNOPT [START_REF] Gill | SNOPT: An SQP algorithm for large-scale constrained optimization[END_REF] , as well as in-house optimization algorithms. 20 DOE algorithms are interfaced from the OpenTURNS 29 and pyDOE libraries. For a given MDO formulation, the driver used to execute the process may be either a DOE or optimization algorithm, and this can be reconfigured easily. The classes hierarchy used to support this feature is shown in Figure 7. This also allows to create trade-off scenarios based on bi-level formulations, where the system level optimization problem is sampled using a DOE algorithm and the sub-scenarios are solved by optimization algorithms.

Visualization of optimization results

MDO results consist in the set of values for the objective function, design variables, coupling variables and constraints, for each iteration of each optimization problem. This data typically represents thousands or millions of scalar values, which cannot be analysed in a single plot. Therefore, a collection of dedicated plots must be used: each one represents a partial view on this data, or projections on adapted spaces. The typical use of these plots is to debug the MDO process, analyse the algorithms performance, or assist decision making in real applications.

In Figure 12, the first graph shows the evolution of the objective value during the optimization. The second graph shows the normalized values of the design variables with the index of the variable in the vector on the Y axis, with respect to the optimization iterations. The third graph shows the distance to the best optimization iteration, in log scale. The last graph shows an approximation of the second order derivatives of the objective function

∂ 2 f (x)
∂x 2 , which is a measure of the sensitivity of the function with respect to the design variables, and of the anisotropy of the problem (differences of curvatures in the design space).

Figure 13 portrays the evolution of the values of the constraints. The components of g 1 , g 2 , g 3 are concatenated as a single vector of size 12. The constraints must be non-positive, that is the plot must be green or white in order to be satisfied (white = active, red = violated). At convergence, only two of them are active : g 2 and g 3 3 .

One may also be interested in getting information on the constraints values while following the objective values iteration per iteration. Figure 14 plots all the constraints components using colors at each iteration while Figure 13 plots the objective function value superimposed with the maximum constraint violation. Therefore the first figure is dedicated to a fine analysis of the constraints values while the latter is used to analyse the best solution in terms of objective and constraints compromise.

In other cases, one may prefer to plot all constraints values per iteration. In Figure 15, the constraints are plotted in separate figures, while colors in the background provide information on the magnitude of the

(

(3) values. Because of the plots scaling, in some cases it can be difficult to identify whether a constraint is slightly violated or not. This is why a vertical black line is plotted at the last iteration before the last sign change of the considered constraint.

(4) (5) (6) (7) g_2 g_3 (1) (2) (3) (4)
All constraints may also be compared at a given iteration in Figure 16 on a radar chart. The scaling allows to compare the values of all constraints simultaneously while having different orders of magnitude. The radar shape makes the plot less sensitive to the number of constraints than classical line plots.

The parallel coordinates portray the design variables history during the scenario execution. In Figure 17, each horizontal coordinate is dedicated to a design variable, normalized by its bounds. A polyline joins all components of a given design vector and is colored by objective function values. This highlights correlations between values of the design variables and the objective function. Figure 17 and Figure 18 show parallel coordinates respectively of the design variables and the constraints, colored by objective function value. The minimum objective value, in dark blue, corresponds to satisfied constraints.

In Figure 19, a Self-Organizing Map (SOM) is generated by using an unsupervised artificial neural network [START_REF] Kohonen | Self-Organizing Maps[END_REF] . A map of size n x × n y is generated, where n x is the number of neurons in the X direction and n y is the number of neurons in the Y direction. The design space (whatever the dimension) is reduced to a 2D representation based on n x n y neurons. Samples are clustered to a neuron when their design variables are close in terms of L2 norm. A neuron is always located at the same place on a map. Each neuron is colored according to the average value for a given criterion. This helps to qualitatively analyze if parts of the design space are good according to some criteria and not for others, and where compromises should be made. SOM provide a qualitative view of the objective and the constraints, and of their relative behaviors. For more details, see the study by Kumano et al. [START_REF] Kumano | Multidisciplinary design optimization of wing shape for a small jet aircraft using kriging model[END_REF] on wing MDO post-processing using SOM. Figure 19 features an SOM example on the Sobieski problem. A DOE of 10 000 samples using the MDF formulation is used to generate the data for the SOM. The process is created using the MDF formulation, that involves an MDA; the objective function and constraints are then called by the DOE algorithm. Figure 20 shows the total derivatives ∂f ∂xi of the objective and constraints with respect to the design variables:

• a large value means that the function is sensitive to the variable,

• a null value means that, at the optimal solution, the function does not depend on the variable,

• a negative value means that the function decreases when the variable is increased.

Figure 20 displays the total derivatives of the MDA with respect to all the design variables at the optimum of the SSBJ test case. x 0 (wing taper ratio) and x 2 (Mach number) appear to be the most important for the gradient of the objective function. Constraints g 1 0 to g 1 4 are very similar, since they all quantify the stress in various sections. g 1 5 and g 1 6 correspond to the lower and upper bounds of the twist, therefore their sensitivities are of opposite signs. g 2 is a function of only x 0 since x 0 is the only variable that influences its gradient.

Figure 21 shows an approximation of the Hessian matrix ∂ 2 f ∂xi∂xj using the Symmetric Rank 1 method (SR1) 32 at the optimum of the SSBJ problem. The color map uses a symmetric logarithmic scale. The Hessian matrix is proportional to the curvature of the function, and non-diagonal terms are the crossinfluence of the variables, which are useful information for understanding the design problem properties.

MDO formulations scalability benchmark

Multiple MDO formulations may be used to solve a given MDO problem. However, their performance vary depending on the addressed problem. The a priori selection of the most appropriate one is not straightforward. To this aim, GEMS supports the scalable MDO problem methodology proposed by Vanaret et al. [START_REF] Vanaret | On the Consequences of the "No Free Lunch" Theorem for Optimization on the Choice of an Appropriate MDO Architecture[END_REF] . This methodology takes a reference set of disciplines and builds surrogate versions of them, using a small DOE. These surrogates have a particular property that allows to change the dimensions of the inputs and output variables as desired. Combined with the automated MDO process generation capabilities of GEMS, one can then assess the relative scalability of the formulations with respect to the dimensions of the problem, without having to solve the real and computationally expensive MDO problem.

Wrapping the disciplines

To interface a simulation program with GEMS, two options are available.

1. Write a generic GEMS wrapper for any software interfaced with a given workflow engine. This option is the best if the workflow engine provides required services, such as HPC access, working directory management, post-processing, or if you want to easily chain multiple programs.

2. Write a specific wrapper for the software by directly inheriting from GEMS' MDODiscipline class in Python (see Section 3). This option is possible if you have no workflow engine, or if you have a light and pure in-memory code that does not require a lot of file processing or multi-machine access.

Both options can be mixed to build MDO processes. JSON schemas 33 are used to describe the inputs, outputs and configuration parameters of most of GEMS components (disciplines, algorithms, visualizations). The JSON schema is a standard for the description of data structures (similar to the XML schema), typically used in web forms. This technology is implemented in most programming languages, which facilitates the development of GUIs to fill data forms and check the data provided by the user. • A generic way to hybridize MDA methods

• A generic way of decomposing MDAs into sub-MDAs

In order to decompose a large coupled problem into multiple sub-coupling problems, the disciplines dependencies are analysed based on the inputs and outputs naming conventions. A directed graph of dependency is created, similar to the N2 matrix. On this graph, two algorithms are applied. First, the Tarjan algorithm 34 determines the strongly connected disciplines. Then, a Kahn's topological sort algorithm [START_REF] Kahn | Topological Sorting of Large Networks[END_REF] determines the steps that may be run in parallel. The sub MDAs may be solved by any of the previously listed algorithms. Coupled derivatives may be computed using the discrete adjoint, as explained in the next section. The overall procedure is summarized in Figure 22.

Coupled derivatives

Gradient-based algorithms require the computation of the total derivatives of the output function f = (f 0 , . . . , f N ) T with respect to the design vector x = (x 0 , . . . , x n ) T . These derivatives are computing while respecting the state equations of the system : R(Y, x) = 0, where R is the whole residual vector of the coupled system and Y is the whole vector of coupling variables. (

GEMS provides several features to compute the derivatives of the objective function and constraints: semi-analytical coupled derivatives, finite differences and complex step [START_REF] Squire | Using Complex Variables to Estimate Derivatives of Real Functions[END_REF][START_REF] Martins | The Complex-Step Derivative Approximation[END_REF] . The semi-analytical derivatives allow significant time savings and higher precision compared to the other options, but it requires that the interfaced disciplines provide the derivatives of their outputs with respect to their inputs, which is not always the case.

When applying finite differences or the complex step method to an MDA (monolithic differentiation), an MDA has to be solved for each perturbed point (x + h j e j ):

df i dx j = f i (x + h j e j ) -f i (x) h j + O(h j ). (2)
If the size of the design vector is large, computing the sensitivity of the output f with respect to the design vector x becomes tedious.

The derivatives computations in GEMS is based on the discrete adjoint method.

• Direct method: dY dx is computed by solving a linear system:

df dx = - ∂f ∂Y • ∂R ∂Y -1 • ∂R ∂x -dY/dx + ∂f ∂x . (3)
• Adjoint method: the computation of the adjoint vector λ is independent of the design variable vector x:

df dx = - ∂f ∂Y • ∂R ∂Y -1 λ T • ∂R ∂x + ∂f ∂x = -λ T • ∂R ∂x + ∂f ∂x . (4)
The choice of the method (direct or adjoint) depends on how the number of design variables n compares to the size of the function (objective and constraints) vector N : if N n, the adjoint method should be used, whereas the direct method should be preferred if n N . Both direct and adjoint methods are implemented in GEMS.

Conclusion

A Python library, GEMS, has been created to automatically build MDO processes, based on a range of MDO formulations, in particular bi-level formulations. The use of abstract concepts makes the process creation fully generic, which on one hand facilitates the construction of nested processes that are key in the implementation of multi-level MDO processes, and on the other hand makes possible to separate the MDO scenario from existing sub-workflows then managed through standard interfaces. Regarding the MDO workflow itself, no workflow engine is needed to implement it. Instead, Python itself is used to build the MDO workflow and dataflow. This provides a high level of flexibility to the process creation, and offers the capability to easily manipulate inputs and outputs, for instance in the purpose to assemble derivatives or to reconfigure the process.

This new software opens interesting perspectives for complex processes creation that embed a high number of disciplines and several levels of sub-processes, either for the needs of future projects, or for higher-fidelity processes.

Figure 1 :

 1 Figure 1: Two steps for building a process in ModelCenter: defining the execution sequence (workflow), and then the data exchanges (dataflow).

Figure 2 :

 2 Figure 2: A data-driven description of a simulation process.

Figure 3 :

 3 Figure 3: A workflow-driven description of a simulation process.

Figure 4 :

 4 Figure 4: The combinatorial issue in MDO process development.

Figure 5 :

 5 Figure 5: Core class diagram of GEMS.

Figure 6 :

 6 Figure 6: GEMS class diagram: optimization process.

Figure 7 :

 7 Figure 7: GEMS class diagram: driver and scenarios.

Figure 8 :

 8 Figure 8: MDA and generic processes class diagram of GEMS.

Figure 9 :Figure 10 :

 910 Figure 9: GEMS execution sequence diagram

Figure 11 :

 11 Figure 11: XDSM of a bilevel MDO process generated from GEMS.

Figure 12

 12 Figure 12: SSBJ MDF optimization history

Figure 13 :

 13 Figure 13: Optimization history: evolution of the inequality constraints in an MDF process.

Figure 14 :

 14 Figure 14: Objective function with respect to the iteration number, with constraint violation information in the background.

Figure 15 :

 15 Figure 15: Constraints values with respect to iteration number, including magnitude of constraint violation information in the background.

Figure 16 :

 16 Figure 16: Constraints values at a given iteration number.

Figure 17 :

 17 Figure 17: Optimization history: parallel coordinates of the design variables in an MDF process for the SSBJ problem.

Figure 18 :

 18 Figure 18: Optimization history: parallel coordinates of the objective function and constraints in an MDF process for the SSBJ problem.

Figure 19 :

 19 Figure 19: A Self-Organizing Map (SOM) in an MDF based process for the SSBJ problem.

Figure 20 :

 20 Figure 20: Coupled derivatives of the objective function and constraints at the optimum point of the SSBJ problem.

Figure 21 :

 21 Figure 21: A Hessian SR1 approximation of the objective function at the optimum point of the SSBJ problem.

4. 7 .

 7 Advanced Multidisciplinary analyses 4.7.1. Available algorithms GEMS has a series of MDA algorithms : • Gauss-Seidel • Jacobi with parallelism • Newton and 10 quasi-newton variants from Scipy.

Figure 22 :

 22 Figure 22: Generation and resolution of an MDA process.

of 26 American Institute of Aeronautics and Astronautics

of 26 American Institute of Aeronautics and Astronautics

Acknowledgments

The authors wish to acknowledge the PIA framework (CGI, ANR) and the members of Technology IRT Saint Exupéry for their support, financial funding and expertise.

x s h a r e d , x 1 , x 2 , x 3 S u b j e c t t o c o n s t r a i n t s : g 1 ( x s h a r e d , x 1 , x 2 , x 3 ) <= 0 g 2 ( x s h a r e d , x 1 , x 2 , x 3 ) <= 0 g 3 ( x s h a r e d , x 1 , x 2 , x 3 ) <= 0 D e s i g n Space : +- ---------+-------------+-------+-------------+------- ---------+-------------+-------+-------------+------- ---------+-------------+---------------------+-------------+-------+ * * * MDO S c e n a r i o r u n t e r m i n a t e d i n 0 : 0 0 : 0 0 . x s h a r e d , x 1 , x 2 , x 3 , y 1 4 , y 3 2 , y 3 1 , y 2 4 , y 3 4 , y 2 3 , y 2 1 , y 1 2 S u b j e c t t o c o n s t r a i n t s : y 3 1 y 3 2 y 3 4 ( x s h a r e d , x 3 , y 2 3 ) = y 3 1 ( x s h a r e d , x 3 , y 2 3 ) -y 3 1 = 0 y 3 2 ( x s h a r e d , x 3 , y 2 3 ) -y 3 2 = 0 y 3 4 ( x s h a r e d , x 3 , y 2 3 ) -y 3 4 = 0 y 2 4 y 2 3 y 2 1 ( x s h a r e d , x 2 , y 3 2 , y 1 2 ) = y 2 4 ( x s h a r e d , x 2 , y 3 2 , y 1 2 ) -y 2 4 = 0 y 2 3 ( x s h a r e d , x 2 , y 3 2 , y 1 2 ) -y 2 3 = 0 y 2 1 ( x s h a r e d , x 2 , y 3 2 , y 1 2 ) -y 2 1 = 0 y 1 2 y 1 4 ( x s h a r e d , x 1 , y 3 1 , y 2 1 ) = y 1 2 ( x s h a r e d , x 1 , y 3 1 , y 2 1 ) -y 1 2 = 0 y 1 4 ( x s h a r e d , x 1 , y 3 1 , y 2 1 ) -y 1 4 = 0 g 1 ( x s h a r e d , x 1 , y 3 1 , y 2 1 ) <= 0 g 2 ( x s h a r e d , x 2 , y 3 2 , y 1 2 ) <= 0 g 3 ( x s h a r e d , x 3 , y 2 3 ) <= 0 D e s i g n Space : +----------+-------------+--------------------+-------------+------- ---------+-------------+--------------------+-------------+------- ---------+-------------+--------------------+-------------+------- ---------+-------------+---------------------+-------------+------- ---------+-------------+---------------------+-------------+------- x s h a r e d S u b j e c t t o c o n s t r a i n t s : g 1 g 2 g 3 ( x s h a r e d ) <= 0 D e s i g n Space : +----------+-------------+-------+-------------+------- ---------+-------------+-------+-------------+------- ---------+-------------+-------+-------------+------- -----+-------------+-------+-------------+------- -----+-------------+-------+-------------+------- -----+-------------+-------+-------------+------- -----+-------------+--------------------+-------------+------- -----+-------------+--------------------+-------------+------- -----+-------------+--------------------+-------------+------- -----+-------------+-------+-------------+------- -----+-------------+-------+-------------+------- -----+-------------+-------+-------------+------- -----+-------------+-------+-------------+------- | f l o a t | +------+-------------+--------------------+-------------+-------+