Siyang Deng
email: siyang.deng@centralelille.fr

Stéphane Brisset

Stéphane Clenet

COMPARATIVE STUDY OF METHODS FOR OPTIMIZATION OF ELECTROMAGNETIC DEVICES WITH UNCERTAINTY

Keywords: Reliability-based design optimization, uncertainty, reliability, safety transformer. 1

This paper compares different probabilistic optimization methods dealing with uncertainties. Reliability-Based Design Optimization is presented as well as various approaches to calculate the probability of failure. They are compared in terms of precision and number of evaluations on mathematical and electromagnetic design problems to highlight the most effective methods.

INTRODUCTION

In most optimization problems, the variables are usually considered as deterministic, i.e. without any variability. This traditional Deterministic Design Optimization (DDO) addresses only the performances but not the reliability and robustness. Since, the manufacturing process, the characteristics of materials, and the dimensions undergo variability, the device performances are in practice not deterministic but uncertain. Therefore, uncertainties related to design parameters are more and more taken into account in all engineering fields.

Various formulations are available in the literature to express optimization problems with uncertainty, which can be mainly divided into Worst-Case Optimization (WCO), Robust Design Optimization (RDO), Reliability-Based Design Optimization (RBDO) and Reliability-Based Robust Design Optimization (RBRDO). WCO is a non-probabilistic approach that is based on minimax problem formulation. For instance, [START_REF] Xiao | Considerations of uncertainty in robust optimisation of electromagnetic devices[END_REF] solves a multi-objective problem that aims to minimize the objective function, its maximum feasible value in a surrounding box, and the greatest component of the objective function derivative in the surrounding box.

The three others are probabilistic approaches that quantify the uncertainty of quantities of interest by probability distribution functions. RDO minimize a weighted sum of the mean value and variance of the objective function subject to deterministic constraints. RBDO minimize the mean value of the objective function subject to constraints on the probability of failure, i.e. constraint violation. Finally, RBRDO [START_REF] Kim | A Comparative Study on Probabilistic Optimization Methods for Electromagnetic Design[END_REF] integrates both last formulations by changing the objective function and constraints at the same time.

A comparative study [START_REF] Ren | Investigation of reliability analysis algorithms for effective reliability-based optimal design of electromagnetic devices[END_REF] of two RBDO Double-Loop Methods (DLM) with Monte Carlo Simulation (MCS) shows the interest of RBDO-DLM compared to MCS that requires a very large sampling to be accurate. However, other RBDO approaches such as Single-Loop (SLM) and Sequential Decoupled Methods (SDM) were not simultaneously investigated and this paper proposes to compare 6 algorithms belonging to the three aforementioned RBDO approaches with MCS in order to highlight the most accurate and the less time consuming. This is performed with a simple mathematical model and the multidisciplinary optimization problem of a safety transformer with uncertainty.

The paper is organized into three parts. Chapter 2 introduces the different categories of RBDO approaches. Two examples are detailed in the chapter 3 and used to compare the different methods. Last chapter is the conclusion.

RELIABILITY-BASED DESIGN OPTIMIZATION

The original formulation without any uncertainty or DDO is expressed as:

(1)

min 𝑑 𝑓(𝑑) 𝑠. 𝑡. 𝑔(𝑑) ≤ 0 𝑑 𝐿 ≤ 𝑑 ≤ 𝑑 𝑈
where 𝑑 is the input design variable, 𝑓(•) and 𝑔(•) are the objective function and the inequality constraint, and 𝑑 𝐿 , 𝑑 𝑈 represent the lower and upper bounds of 𝑑, respectively.

As the variability of the design variables is taken into account, the original deterministic input parameter 𝑑 should be replaced by a random input parameter 𝑋, which follows the normal law in this paper for simplicity. The mean value of 𝑋 is denoted 𝑑 and is the unknown of the new design variables, the standard deviation is denoted 𝜎 and is considered constant which means that the variability of the input parameter X doesn't depend on its mean.

The method RBDO aims to find the optimal design with the allocation of a target reliability level. The probability of failure is close to the number of points in a sampling around the mean value that fall within the reliable domain. This means that RBDO attempt to find the optimal design with the reliability that the probability of failure must be smaller than a given target value. The formulation of RBDO is as follows:

(2)

min 𝑑 𝑓(𝑑) 𝑠. 𝑡. 𝑃 𝑓 (𝑔(𝑥) > 0) ≤ 𝑃 𝑡 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎
where 𝑃 𝑓 is the probability of failure, 𝛽 𝑡 is the reliability index, and 𝑃 𝑡 = 𝛷(-𝛽 𝑡 ) the target value for the probability of failure. The RBDO uses probabilistic constraints to make sure that the design variables satisfy a desired reliability level while minimizing the function objective.

Probability of failure

The probability is the likelihood of an event, estimated by a real number between 0 and 1. The Probability Density Function (PDF) and the Cumulative Distribution Function (CDF) define the occurrence of stochastic quantities inherently uncertain. The statistical description of a random variable 𝑋 given by the CDF Fx or PDF fX is expressed as follows:

(3)

𝑃[𝑋 ≤ 𝑥] = 𝐹 𝑋 (𝑥) = ∫ 𝑓 𝑋 (𝜏)𝑑𝜏 𝑥 -∞
where 𝑃 is the probability of occurrence of an event.

In the field of electromagnetic device manufacturing, the system ability to satisfy consumer's demand or operating constraints is important. The reliability means that designers should reduce the probability of failure as much as possible. The determination of the reliability of the system is based on the limit state function. Each constraint 𝑔(𝑋) ≤ 0 can separate the domain of 𝑋 into three parts: the limit state curve is 𝑔(𝑋) = 0, the domain where 𝑔(𝑋) > 0 is the failure domain, on the contrary the security domain represents the area 𝑔(𝑋) < 0. where 𝐷 𝑓 is the failure domain. Because the computational burden is heavy with numerous random parameters, direct integration is almost impossible and thus Monte-Carlo Simulation (MCS) or other techniques such as First-Order Reliability Method (FORM) [START_REF] Liu | Optimization algorithms for structural reliability[END_REF], [START_REF] Hasofer | Exact and invariant second-moment code format[END_REF] is often used to calculate an approximation of 𝑃 𝑓 .

FORM and inverse FORM are based on an isoprobalistic transformation to have a normalized vector of statistically independent random variables 𝑈 instead of the initial input parameter 𝑋.

For the Gaussian vector 𝑋 in this paper, the transformation 𝑇 is as follows: [START_REF] Hasofer | Exact and invariant second-moment code format[END_REF] 𝑇: 𝑈 = (𝑋 -𝑑) 𝜎 ⁄ Figure 2. The first order reliability method [START_REF] Lopez | Reliability-based design optimization strategies based on FORM: a review[END_REF].

Then the limit state function changes from 𝑔(𝑑) = 0 to 𝐺 𝑈 (𝑢) = 0 and the CDF from 𝐹 𝑋 to 𝐹 𝑈 . 𝐺 𝑈 (𝑢) is defined as the performance function, and the FORM or inverse FORM method uses a linear approximation to replace the real performance function at the Most Probable Point of failure (MPPF) in U-space. The MPPF 𝑢 * is the one that minimize the distance between the origin 𝑂 and 𝐺 𝑈 (𝑢) = 0. After 𝑢 * is found, the limit-state function is replaced by a tangent hyperplane crossing 𝑢 * . So the probability of failure is calculated by: ( 6)

𝑃 𝑓 ≈ ∫ 𝑓 𝑈 (𝑢)𝑑𝑢 = 𝛷(-𝛽) 𝐺 ̃𝑈(𝑢)>0
where 𝐺 ̃𝑈(𝑢) = 0 is the hyperplane that approximates the limit-state function 𝐺 𝑈 (𝑢) = 0, 𝛷 is the standard Gaussian cumulative distribution function, and 𝛽 = ‖𝑢‖ is the reliability index which is equal the norm of 𝑢.

The formulation ( 2) is the basis of all RBDO methods and can be solved by different approaches which are usually separated in three main categories: double-loop, single-loop, and sequential decoupled methods. The following sections introduce the principle of these different categories and present some approaches for each of them.

Double-loop methods

Double-loop methods use two loops to solve RBDO problems: the inner loop aims to analyze reliability of the chosen configuration and to calculate the probability of failure using FORM or inverse FORM; the outer loop seeks the mean values of input designs that minimize the objective function and constrain the probability of failure computed by inner loop.

There are several approaches for double-loop methods, the most popular are Reliability Index Approach (RIA) and Performance Measure Approach (PMA) [START_REF] Ren | Investigation of reliability analysis algorithms for effective reliability-based optimal design of electromagnetic devices[END_REF], [START_REF] Tu | A new study on reliability-based design optimization[END_REF].

Reliability Index Approach (RIA)

RIA uses the FORM to calculate the reliability index in the inner loop:

(7) 𝛽 = min 𝑢 ‖𝑢‖ 𝑠. 𝑡. 𝐺 𝑈 (𝑢) = 0
The outer loop of RIA, any constrained non-linear algorithm like SQP or others can be chosen to minimize the objective function 𝑓 with the index 𝛽 given by the FORM:

(8) min 𝑑 𝑓(𝑑) 𝑠. 𝑡. 𝑔 ≤ 0 𝛽 ≥ 𝛽 𝑡 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎
where 𝑔 < 0 is used to restrict the because the definition in equation ( 6) only validate if the origin is located in the security domain.

Performance Measure Approach (PMA)

In former formulation of RBDO, optimization is carried out with the limitation of the reliability index that must be greater than or be equal to a target value. The calculation of this index leads to the search for the MPPF. In contrast, for PMA formulation, optimization is formulated with the limitation of maximum performance with a given value of reliability index. The search of this maximum performance is to maximize the function 𝐺 𝑈 with the limitation that reliability index must achieve the target value. This approach is considered as the inverse of the FORM approximation [START_REF] Wu | Advanced probabilistic structural analysis method for implicit performance functions[END_REF], [START_REF] Youn | Hybrid analysis method for reliability-based design optimization[END_REF].

So the outer loop becomes: where 𝑢 * is the Maximum Performance Target Point (MPTP) that corresponds to the target index 𝛽 𝑡 .

Single-loop methods

For the so-called single-loop methods, the main point is that the inner loop is replaced by an approximation to avoid the iterative evaluations for reliability analysis in order to accelerate the convergence to the optimum.

Approximate Moment Approach (AMA)

AMA is based on statistical moments. The first order Taylor expansion is used to calculate the mean 𝜇 𝑔 𝑖 and the standard deviation 𝜎 𝑔 𝑖 for all constraint functions 𝑔 𝑖 (𝑑) by using the following expressions [START_REF] Putko | Approach for input uncertainty propagation and robust design in CFD using sensitivity derivatives[END_REF]:

(11) { 𝜇 𝑔 𝑖 = 𝑔 𝑖 (𝑑) 𝜎 𝑔 𝑖 2 ≈ (𝛻 𝑔 𝑖 𝑇 ) 2 • 𝜎 2
where 𝜎 is the standard deviation of the input variables 𝑋 and 𝛻 is the gradient operator. With these expressions, an approximation 𝑔 𝑖 (𝑑) + 𝛽 𝑡 𝜎 𝑔 𝑖 is used to replace the original constraints 𝑔 𝑖 (𝑑′) as shown in Figure 3. In this way, the margin of security is shifted in order to keep the reliability as desired. Figure 3 shows the basic idea of AMA approach, 𝑔(𝑑) = 0 is the deterministic limit-state function. As the optimum point 𝑑 * is on the limit-state curve, the probability of failure is about 50%. The situation which could cause the failure around this point with a given probability is shown as shaded area on the left side. This approach is based on a local linear approximation of the constraint functions at the mean value of the design parameters and as the probabilistic distribution of the performances depends on two moments that are also approximated, the reliability assessment could produce significant numerical error [START_REF] Youn | Selecting probabilistic approaches for reliability-based design optimization[END_REF].

Single Loop Approach (SLA)

SLA uses the same approximation as AMA for calculating the moments and the position of MPPF (marked as 𝑑 ′ * in Figure 3 and 𝑥 𝑝 in Figure 4). The difference between these two methods is that SLA evaluates the constraints 𝑔 at the approximate MPPF unlike AMA that uses a first-order approximation of the constraint function around the mean value 𝑑 directly. The formulation is given as [START_REF] Liang | A single-loop method for reliability-based design optimization[END_REF]: where 𝛼 𝑖 is the normalized gradient of the 𝑖th constraint, 𝑥 𝑝 𝑖 is the approximation of the MPPF for the constraint 𝑔 𝑖 , and ∘ is the Hadamard operator (element-wise) product.

Figure 4. Principle of the SLA approach

The principle of this approach is similar to AMA. As shown in Figure 4, it does not search for the MPPF of each constraint by using an inner loop but approximate its position.

Sequential decoupled methods

Sequential decoupled methods aim to change the initial problem into a series of optimization cycles. The cycles are sequential, each individual optimization is deterministic and uses the optimum given by the former optimization as an initial point. At the first iteration, the algorithm searches a deterministic solution without considering uncertainty and then compute the reliability index of this solution to deduce a shift in order to achieve a given probability of failure. The next iterations refine the shift.

Sequential Optimization and Reliability Assessment (SORA)

SORA employs a series of cycles of deterministic optimizations and reliability assessments. In each cycle, optimization and reliability assessment are decoupled from each other, the reliability assessment is only conducted after the deterministic optimization to verify constraint feasibility under uncertainty. The main point of this method is to shift the boundaries of violated constraints to the feasible direction based on the reliability information obtained in the former cycle [START_REF] Du | Sequential optimization and reliability assessment method for efficient probabilistic design)[END_REF]. The updated point is used in the next cycle of the deterministic optimization. This cycle is repeated until the fulfillment of the convergence criteria.

The process of SORA method is presented as follows. First of all, an initial shift 𝑠 0 = 0 allows finding the solution 𝑑 * 0 of the deterministic problem. At each iteration, the optimization problem is defined by: where 𝑘 indicates the iteration number. The optimal value is set as 𝑑 * 𝑘 . After each optimization, some of the constraints may become active. For an active constraint, the optimal point 𝑑 * 𝑘 is on the limit-state curve. When considering the randomness of 𝑋, the actual probability of failure is about 0.5, so a reliability assessment is implemented at the deterministic optimum solution to locate the MPTP 𝑥 𝑘 that corresponds to the desired probability of failure. To ensure the MPTP onto the deterministic boundary, a shifting vector 𝑠 𝑘+1 is derived: [START_REF] Tran | A Benchmark for Multi-Objective, Multi-Level and Combinatorial Optimizations of a Safety Isolating Transformer[END_REF] 𝑠 𝑘+1 = 𝑑 * 𝑘 -𝑥 𝑘 Therefore, when establishing the equivalent deterministic optimization model in the next cycle, the constraints is modified to shift the MPTP at least onto the deterministic boundary.

Sequential Approximate Programming (SAP)

SAP is another approach based on Taylor expansion at the first order. The original optimization problem is decomposed into a sequence of sub-optimization problems. Each sub-optimization consists of an approximate objective function subjected to a set of approximate constraint functions [START_REF] Yi | A sequential approximate programming strategy for performance-measure-based probabilistic structural design optimization[END_REF]. The details are presented as follows:

For considering the PMA formed optimization problems, the expressions are as equation ( 8), following the idea of SAP, a sequential approximate formulation is constructed as: where 𝐺 ̂𝑝(𝑑 𝑘-1 ) = 𝐺 𝑢 (𝑢 𝑘-1 ) is the approximate probabilistic performance measure and 𝛻 𝑑 𝐺 ̂𝑝(𝑑 𝑘-1 ) is its gradient. To avoid another optimization to calculate 𝑢 𝑘 , it is updated by the function below:

(19) 𝑢 𝑘 = -𝛽 𝑡 𝛻 𝑈 𝐺 𝑈 (𝑢 𝑘-1 )
‖𝛻 𝑈 𝐺 𝑈 (𝑢 𝑘-1 )‖

where 𝑢 0 = 0 is usually chosen as the initial estimation.

With the first order approximation of 𝐺 𝑝 (𝑑 𝑘 ), new formulation is obtained:

(20)

min 𝑑 𝑓 𝑘 (𝑑) 𝑠. 𝑡. 𝐺 ̂𝑝(𝑑 𝑘-1 ) + (𝛻 𝑑 𝐺 ̂𝑝(𝑑 𝑘-1 )) 𝑇 (𝑑 -𝑑 𝑘-1 ) ≤ 0 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎
This formulation of SAP is similar to the outer loop of PMA. The same principle is used to convert a RIA formulation into a sequential approximate programming [START_REF] Cheng | A sequential approximate programming strategy for reliability-based structural optimization[END_REF].

CASE STUDIES

In this chapter, two examples are tested to compare the performances of all approaches.

Numerical example

To assess the efficiency of these methods, the simple two-variable problem in [START_REF] Kim | A Single-Loop Strategy for Efficient Reliability-Based Electromagnetic Design Optimization[END_REF] is analyzed. The optimization problem of this numerical example is: In order to understand the implications of the choice of different formulations and also to quantify the impact of the choice of approaches and algorithms on the accuracy and the number of evaluations, the results of the 7 mentioned methods are given in Table 1. The different optima are compared with each other in Figure 5, where the dotted curves present contours of 𝑓 and constraint boundaries 𝑔 𝑖 = 0 are depicted as solid lines. The target probability of failure 𝑃 𝑡 for RBDO is 2.28%, 𝑃 𝑓 is computed by Monte-Carlo simulation with a sampling of 10 6 realizations. For the simplicity and without loss of generality, all the uncertainties are modeled with the normal law. The column named convergence in Table 1 means that we run the algorithms 100 times with different initial points to see how many times it converges to the same optimum. From Table 1 and Figure 5, we can see that most of the methods can find a result that satisfies the constraints. DDO leads to a solution with the lowest objective function mean value, the highest standard-deviation of objective function and a probability of failure around 50%. On the contrary, the probabilities of failure of RBDO methods are close to the target probability except for the single-loop method AMA because of the approximation used to increase the speed of convergence by sacrificing the accuracy. Among RBDO methods, single-loop strategies are fast but inaccurate, double-loop and sequential decoupled methods lead to the same results but sequential decoupled methods are greatly faster.

Example of safety transformer optimization

The example of a safety isolating transformer [START_REF] Tran | A Benchmark for Multi-Objective, Multi-Level and Combinatorial Optimizations of a Safety Isolating Transformer[END_REF] is also investigated. This is a single-phase transformer with grain-oriented E-I laminations. The primary and secondary windings are both wound around the frame surrounding the central core. The model inputs consist in 7 random design parameters: three parameters 𝑎, 𝑏, 𝑐 for the shape of the lamination, one for the frame 𝑑, two for the section of conductors 𝑆 1 , 𝑆 2 and the last one for the number of primary turns 𝑛 1 . These design variables are shown in Figure 6. The range of input variables are:

(22) 3𝑚𝑚 ≤ 𝑎 ≤ 30𝑚𝑚 14𝑚𝑚 ≤ 𝑏 ≤ 95𝑚𝑚 6𝑚𝑚 ≤ 𝑐 ≤ 40𝑚𝑚 10𝑚𝑚 ≤ 𝑑 ≤ 80𝑚𝑚 200 ≤ 𝑛 1 ≤ 1200 0.15𝑚𝑚 2 ≤ 𝑆 1 ≤ 19𝑚𝑚 2 0.15𝑚𝑚 2 ≤ 𝑆 2 ≤ 19𝑚𝑚 2
There are 8 inequality constraints in this problem. The copper and iron temperatures 𝑇 𝑐𝑜𝑛𝑑 and 𝑇 𝑖𝑟𝑜𝑛 should be less than given temperatures 120℃ and 100℃ respectively. Both the magnetizing current should be less than 10%. The filling factors of the primary coil 𝑓 1 and secondary coil 𝑓 2 should be lower than 1, and the efficiency 𝜂 should be larger than 0.8. At last, the residue must be less than 10 -6 .

The constraint functions are: The aim is to minimize the total mass 𝑀 𝑡𝑜𝑡 with a target probability 𝑃 𝑡 = 0.13%, so that the target reliability index is 𝛽 𝑡 = 3.

Table 2 shows the results of some methods. For each constraint, there are two probabilities of failures, upper ones are calculated by the methods themselves and the below ones are calculated by the Monte-Carlo simulation for comparison. For the reason that the 8 th constraint is a residue, we consider it as a deterministic one, so there are only 7 probabilities of constraints presented in the table 2. Note that PMA can find an optimum that satisfies all the constraints even if the number of evaluations of the model is high. For SLA, it has a smaller number of evaluations but there is one constraint violated. The reason is that SLA sacrifices the accuracy in order to reduce the number of evaluations, leading to a coarse computation of the probability of failure. The convergence rate of SORA is not as good as the other two but it has the smallest number of evaluations among them. It can be seen that this number is nearly 15 times less than PMA and 3 times less than SLA. Unfortunately, the rate of convergence for SORA is 4 times lower. This means that a multi-start process is required and the number of evaluations will increase consequently. We obtain almost the same conclusions as with the numerical example: the single-loop method SLA is the most inaccurate method and doubleloop method PMA has the highest convergence rate. For this more complicated example, the number of evaluations of sequential decoupled method SORA is even smaller than of SLA, so that the more complex the model is, the more effective SORA may be. Other methods fail to find a solution, probably because this problem is hardconstrained and the solution of DDO is on the limit-state of four constraints. So, it also indicates that not all the aforementioned approaches can handle complicated models.

CONCLUSION

RBDO methods change the initial constraints into probabilistic ones and use different approaches to approximate the limit state function or the most probable point of failure to calculate the probability of failure or the reliability index.

The mathematical example shows that the 6 RBDO approaches have almost the same results except AMA that is less accurate. The optimization of the safety transformer highlights that not all the methods can converge to the global solution. PMA, SLA, and SORA appear to be more stable. Considering both numerical examples, SORA is the most effective method among all RBDO approaches.

Figure 1 .

 1 Figure 1. The failure domain, security domain, and limit-state curve. The probability of failure is the probability of the event 𝑔(𝑋) > 0. It is calculated with the integral: (4) 𝑃 𝑓 = 𝑃[𝑔(𝑋) > 0] = 𝑃[𝑋 ∈ 𝐷 𝑓 ] = ∫ 𝑓 𝑋 (𝑥)𝑑𝑥 𝐷 𝑓

  . 𝐺 𝑝 (𝑢 * ) ≤ 0𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎where 𝐺 𝑝 (𝑢 * ) is the maximal performance measurement. The purpose of the inner loop is to find 𝑢 * in U-space such as:[START_REF] Youn | Selecting probabilistic approaches for reliability-based design optimization[END_REF] 𝐺 𝑝 (𝑢 * ) = max 𝑢 𝐺 𝑈 (𝑢) 𝑠. 𝑡. ‖𝑢‖ = 𝛽 𝑡

Figure 3 .

 3 Figure 3. Principle of the AMA approach

  . 𝑔 𝑖 (𝑥 𝑝 𝑖 ) ≤ 0 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎 with (14) 𝑥 𝑝 𝑖 = 𝑑 -𝛽 𝑡 𝛼 𝑖 ∘ 𝜎 𝛼 𝑖 = 𝛻𝐺 𝑈 𝑖 (𝑑 ) ‖𝛻𝐺 𝑈 𝑖 (𝑑 )‖ = 𝜎∘𝛻𝑔 𝑖 (𝑑 ) ‖𝜎∘𝛻𝑔 𝑖 (𝑑 )‖

  . 𝑔(𝑑 𝑘 -𝑠 𝑘 ) ≤ 0 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎

  . 𝐺 𝑝 (𝑑 𝑘 ) ≤ 0 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎The approximate function 𝐺 𝑝 (𝑑 𝑘 ) is built with a first order Taylor expansion with respect to the design variables at the current point:(18)𝐺 𝑝 (𝑑 𝑘 ) ≈ 𝐺 ̂𝑝(𝑑 𝑘-1 ) + (𝛻 𝑑 𝐺 ̂𝑝(𝑑 𝑘-1 )) 𝑇 (𝑑 𝑘 -𝑑 𝑘-1 )

Figure 6 .

 6 Figure 6. 3D geometric simulation of the transformer

Table 1 . Results for different optimizations

 1

	Methods	𝑑 1	𝑑 2	𝜇 𝑓	𝜎 𝑓	Max of 𝑃 𝑓 (%)	Convergence (100 times)	Function evaluations of 𝑓 and 𝑔
	DDO	2.4398	0.8400	-2.6266	0.2072	50	20%	84
	RIA	2.2513	1.9691	-1.9945	0.1792	2.28	29%	4539
	PMA	2.2513	1.9691	-1.9945	0.1792	2.28	82%	3183
	AMA	2.3524	2.1568	-1.8712	0.1715	1.87	46%	132
	SLA	2.2512	1.9677	-1.9953	0.1793	2.29	32%	165
	SORA	2.2513	1.9691	-1.9945	0.1792	2.28	45%	531
	SAP	2.2513	1.9691	-1.9945	0.1792	2.28	37%	181

Figure 5. The results of different methods plotted in search space.

Table 2 . Results on the safety transformer optimization problem

 2

		Transformer		DDO	PMA	SLA	SORA
		𝑎 (𝑚𝑚)		13.004	13.078	13.045	13.078
		𝑏 (𝑚𝑚)		50.1	52.274	51.613	52.278
		𝑐 (𝑚𝑚)		16.537	16.758	16.837	16.757
		𝑑 (𝑚𝑚)		43.05	42.182	42.603	42.180
		𝑛 1 𝑆 1 (𝑚𝑚 2 ) 𝑆 2 (𝑚𝑚 2 )		639.76 0.3238 2.9026	658.36 0.3256 2.9296	653.56 0.3254 2.9273	658.39 0.3256 2.9296
		𝜇 𝑀 𝑡𝑜𝑡		2.3112	2.3546	2.3534	2.3546
		𝜎 𝑀 𝑡𝑜𝑡		0.0079	0.0081	0.0081	0.0081
		Convergence		0.12	0.84	0.72	0.23
	𝑃 𝑓 1 (𝑇 𝑐𝑜𝑛𝑑 > 120℃) (%)	by MC by itself	0 0	0 0	0 0	0 0
	𝑃 𝑓 2 (𝑇 𝑖𝑟𝑜𝑛 > 100℃) (%)	by MC by itself	50 50.658	0.13 0.1337	0.13 0.3754	0.13 0.1389
	𝑃 𝑓 3 (∆𝑉 2 𝑉 20 ⁄	> 0.1) (%)	by MC by itself	0 0	0 0	0 0	0 0
	𝑃 𝑓 4 (𝐼 10 𝐼 1 ⁄ > 0.1) (%)	by MC by itself	50 50.203	0.13 0.1356	0.13 0.1357	0.13 0.1370
	𝑃 𝑓 5 (𝑓 1 > 1) (%)	by MC by itself	50 49.987	0.13 0.1323	0.13 0.1362	0.13 0.1388
	𝑃 𝑓 6 (𝑓 2 > 1)(%)	by MC by itself	50 50.021	0.13 0.1415	0.13 0.1368	0.13 0.1382
			by MC	50	0.13	0.13	0.13
			by itself	49.757	0.1324	0.1376	0.1390
	𝑃 𝑓 7 (0.8 -𝜂)(%)
	Model evaluation		242	30607	6048	2171


Laura Balea

Gilles Dusserre
email: gilles.dusserre@mines-albi.fr

Gérard Bernhart

Mechanical behaviour of plain-knit reinforced injected composites: Effect of inlay yarns and fibre type

Keywords: A, Fabrics/textile B, Mechanical properties

This paper investigates the effect of inlay yarns and fibre type (E-glass, basalt, carbon) to improve the mechanical behaviour of plain knit reinforced composites (epoxy matrix). The tensile behaviour of the dry reinforcement was investigated in the wale and course-directions and has shown that the coursewise deformation is drastically reduced whereas the strength is strongly increased when inlay yarns are included. Fibre type has a weak effect on deformation and strength. For the composite material, processed by Liquid Composite Moulding, the carbon fibre reinforcement gives always the best results whatever the testing direction. Inlay yarns decrease the 0°mechanical properties, whereas 45°and 90°are increased. Moreover a quasi-isotropic behaviour is obtained when using two inlay yarns, irrespective of the fibre types.

Introduction

Fibre reinforced composite materials are increasingly used in industrial applications requiring high specific mechanical properties. They cover a large range of mechanical and physical properties due to the numerous parameters involved in their manufacturing. Amongst them, the fibre arrangement in the matrix has an important effect on the mechanical properties of the composite. Up to now, woven and unidirectional reinforcements are mainly used in industrial applications. The low curvature of the fibres confers high performance to the composite. However, issues related to fibrous architecture remain, especially for composite forming. Complex shape parts manufacturing requires cutting and assembly operations, resulting in material wastage and high production costs. This leads also to a local weakening in assembly areas, mainly related to the reinforcement discontinuity. Therefore new reinforcement architectures were designed, and, amongst them, knitted fabrics meet the technical requirements for complex shapes forming. This technology enables the development of netshape fabrics, including holed preforms. Knits are obtained by intermeshing loops of fibre yarns resulting in a high deformation ability compared with woven fabrics [START_REF] Luo | Biaxial tension and ultimate deformation of knitted fabric reinforcements[END_REF][START_REF] Gommers | Analysis of knitted fabric reinforced composites: Part 1. Fibre orientation distribution[END_REF] thanks to the possibility to change the shape of the loops [START_REF] Araújo | Modelling and simulation of the mechanical behaviour of weft knitted fabrics for technical applications -Part 1: General considerations and experimental analyses[END_REF][START_REF] Hagège | Finite element analyses of knitted composite reinforcement at large strain[END_REF]. This deformation ability is especially suitable for stretch forming and deep drawing, since it allows to distribute homogeneously the textile deformation, avoiding local material heterogeneities in a single part [START_REF] Lim | Optimization of the formability of knitted fabric composite sheet by means of combined deep drawing and stretch forming[END_REF][START_REF] Lim | Axisymmetric sheet forming of knitted fabric composite by combined stretch forming and deep drawing[END_REF]. A low strain energy is required allowing wrinkles to appear lately in comparison to woven fabrics [START_REF] Rozant | Drapability of dry textile fabrics for stampable thermoplastic preforms[END_REF][START_REF] Rozant | Warp-knit laminates for stampable sandwich preforms[END_REF]. Excellent properties for composite moulding are obtained, especially for stacking.

Knit reinforced composites exhibit also interesting impact behaviour properties. Stitches of adjacent layers overlap each other, reducing inter-laminar resin-rich areas and crack propagation. The inter-laminar fracture toughness is then increased comparing to woven and UD composites [START_REF] Khondker | Impact and compression after-impact performance of weft-knitted glass textile composites[END_REF][START_REF] Mouritz | Mode I interlaminar fracture toughness properties of advanced textile fibreglass composites[END_REF][START_REF] Falconnet | Fracture toughness of weft-knitted fabric composites[END_REF][START_REF] Pandita | Impact properties of weft knitted fabric reinforced composites[END_REF]. Indeed the delamination of knit composites involves both mode I and mode III crack propagation [START_REF] Pandita | Impact properties of weft knitted fabric reinforced composites[END_REF]. Compression-After-Impact tests highlight the high damage tolerance of knitted composites: they exhibit a relative compressive strength retention higher than woven composites [START_REF] Khondker | Impact and compression after-impact performance of weft-knitted glass textile composites[END_REF].

However, knitted composites exhibit low in-plane mechanical properties, close to short fibres composites, because of the curvature of the fibres [START_REF] Verpoest | The potential of knitted fabrics as a reinforcement for composites[END_REF][START_REF] Huang | Micromechanical modeling approaches for the stiffness and strength of knitted fabric composites: a review and comparative study[END_REF]. The compressive properties are isotropic and depend strongly on the matrix properties [START_REF] Leong | An investigation of the mechanical performance of weft-knit milano-rib glass/epoxy composites[END_REF]. On the contrary, the tensile properties are anisotropic and related to the reinforcement. Up to now, most studies were focused on basic knitted structures such as plain, 1 Â 1 rib and Milano made of glass fibres. These investigations highlighted that the composite properties depend on knit architecture [START_REF] Khondker | An investigation of the structureproperty relationship of knitted composites[END_REF][START_REF] Abounaim | Development of flat knitted spacer fabrics for composites using hybrid yarns and investigation of twodimensional mechanical properties[END_REF] and on knit geometrical parameters such as stitch size, stitch density and yarn length [START_REF] Khondker | An investigation of the structureproperty relationship of knitted composites[END_REF][START_REF] Anwar | The effect of architecture on the mechanical properties of knitted composites[END_REF]. An initial fabric deformation before composite processing leads to an improvement of the properties in the stretched direction and to a loss of properties in the transverse direction [START_REF] Leong | The effects of deforming knitted glass fabrics on the basic composite mechanical properties[END_REF]. However biaxial stretching before composite forming has a negligible effect [START_REF] Khondker | Effects of biaxial deformation of the knitted glass preform on the in-plane mechanical properties of the composite[END_REF]. It is then demonstrated that an improvement of the tensile properties would require a modification of the fibre distribution in the material. No studies on the effect of the fibre type on the knitted composite behaviour are known.

In this paper, the effects of knit construction and fibre type on the composite tensile behaviour are investigated. Previous works [START_REF] Balea | Mechanical behaviour of technical knitted fabrics for composite applications[END_REF][START_REF] Dusserre | Influence of knitted reinforcement on the composite mechanical behaviour[END_REF] were focused on the two basic knitted reinforcements, plain and 1 Â 1 rib. This latter is made up of front and back stitches, leading to a three-dimensional structure and a higher deformability is then achieved than with plain-knit, only made up of front stitches. Plain structure confers the best tensile properties to the composite, especially in the wale-direction, due to a less isotropic fibre orientation distribution, and has been chosen as the starting point of the optimisation aiming to improve the composite mechanical properties [START_REF] Balea | Comportement des matériaux composites à renforts tricotés élaborés par injection de résine[END_REF][START_REF] Balea | Mechanical behaviour of technical knitted fabrics for composite applications[END_REF][START_REF] Balea | Weft-knitted composites: Relationships between fabric architecture, fibre type and mechanical properties of the infused composite[END_REF]. Plain-knit composites are anisotropic with higher properties in the wale-direction. These reinforcements have higher deformation ability in the coursedirection. A correlation between fabric deformation ability and composite properties was carried out. The optimisation of the plain-knit architecture requires to reinforce the course-direction and to find a compromise with the deformation ability (reinforcing means adding straight yarns, what obviously decreases the extensibility). The knitting process allows to add non-knitted yarns in the course-direction into the fabric, thanks to inlay yarns obtained by two techniques: float loops [START_REF] Leong | The potential of knitting for engineering composites -a review[END_REF] or tuck loops [START_REF] Leong | The potential of knitting for engineering composites -a review[END_REF][START_REF] Khondker | Initial fracture of the welt weftknitted textile composites[END_REF]. In this paper, the cross effects of inlay yarns, their number, and the fibre type on knitted fabric deformability and composite mechanical behaviour are investigated.

Materials and process

Fibres and resin

E-glass, basalt and carbon fibres are studied to highlight the effect of the fibre type on the behaviour of the fabric and the composite. Glass and carbon fibres are the most widely used in the composite manufacturing industry. The former are mostly dedicated to mass production materials due to its excellent price-quality ratio and the latter provide high cost structural materials, with excellent mechanical properties. Basalt fibres, considered as emergent, have interesting physicochemical properties and slightly better mechanical properties than glass fibres.

The yarn diameter is a key parameter of a knit, it controls the jamming of the loops [START_REF] Balea | Mechanical behaviour of technical knitted fabrics for composite applications[END_REF]. In order to produce equivalent knitted fabrics, roving bobbins with a similar yarn diameter were chosen. This diameter is calculated from diameter and number of filaments by considering a hexagonal fibre packing. Fibres and roving properties are given in Table 1. Glass and basalt fibres have similar properties and both rovings have the same linear density. However a smaller filament diameter confers to basalt yarns a lower bending stiffness. Carbon fibres have thinner filaments and exhibit a higher Young's modulus and a lower density, that confer to the yarns lower bending stiffness and linear density.

The resin used is a two components epoxy system (SICOMIN SR 1710 resin with SD 8822 hardener) designed for resin transfer moulding processes. This system is characterized by a very low viscosity and high mechanical properties, conferring a high inter-laminar shear strength to the composite (Table 1).

Knitted fabrics

Knitting consists in intermeshing loops of yarns to produce fabrics from a single or from several yarns. Weft-knitting refers to the knitting direction: the stitches are knitted successively in the weftdirection, also called course-direction. The loops are interconnected in the wale direction called also columns. Plain-knit, only made of front stitches, is the basic architecture (Fig. 1) in which inlay yarns were added to reinforce the course-direction. Every four plain stitches, inlay yarns are linked to the knit by tuck stitches (Fig. 2). Between two consecutive tuck stitches, the inlay yarn forms a float loop. Basic plain-knit (P0) and plain-knits with one (P1), two (P2) and three (P3) inlay yarns per course were manufactured and are presented in Figs. 1 and2.

All knits under study were produced from each fibre type with the same knitting parameters in order to get equivalent stitch sizes, given by the geometrical parameters defined in Fig. 1. W, C and L Y are respectively the average wale width, course height and yarn length per stitch. W and C were calculated from fabric dimensions and number of stitches, whereas L Y is calculated from fabric weight and yarn linear density. The mean values calculated for all knits under study are W = 4.2 ± 0.4 mm, C = 2.5 ± 0.1 mm and L Y = 17.5 ± 3.9 mm. The important dispersion of these parameters (5-20%) is related to the properties of the different fibres. The stitch size depends on knitting parameters but also on the loop 2. The thickness is measured under a 5 kPa pressure following the Bueno et al. [START_REF] Bueno | Influence of properties at micro and mesoscopic levels on macroscopic level for weft knitted fabrics[END_REF] conditions. For similar yarn diameters, the thickness mainly depends on the fabric texture and less on the fibre type. Indeed the variations are due to the fibre rearrangement, facilitated by a lower diameter and bending stiffness of carbon fibres. The areal density is similar for glass and basalt knits, and lower for carbon fabrics, in agreement with their yarn linear density.

Composite laminates

The reinforcements aforementioned were processed by VARTM in an aluminium rigid closed mould coated with a semi-permanent release agent (Loctite 44-NC Frekote). A vacuum pressure of 1 mbar was applied into the mould by vents located at each corner and the resin was injected through an inlet located at the centre of the mould (Fig. 3). The resin injection was performed at room temperature and the laminates were cured at 60 °C for 16 h on a heating plate. The composite laminates characteristics are given in Table 3. P0 composites were manufactured with 3 reinforcing layers to reach a fibre volume fraction close to 30%. P1 composites were manufactured with two layers providing the same fibre content as P0. P2 and P3 composites contain also two layers as P1, but exhibit a higher fibre volume fraction. The porosity ratio varies between 2.6% and 10.1%. It is lower in P0 laminates, but no correlation can be made with the number of inlays. High ratios are mainly related to air bubbles entrapment during resin mixture. A careful selection of the specimen was carried out, so that the materials tested have the less voids as possible. dard dedicated to uniaxial tensile testing of woven fabrics, since no standards concerning knitted fabrics are available. The fabrics were directly knitted at the testing sample size (300 Â 100 mm 2 ) since the knit cannot be cut without unravelling. A few courses of an elastomeric fibre were knitted on the long edges of the course-direction specimens to avoid unravelling while limiting its intrusiveness. The tests were carried out on an INSTRON 5800R electromechanical testing machine fitted out with pneumatic grips at the constant crosshead displacement rate of 50 mm min À1 . An acceptable dispersion of the experimental data made it possible to select a representative load versus strain curve for the discussion, among the three available for each test condition. The results are given in load per wale or course versus true strain (Fig. 4).

Mechanical behaviour of the knitted fabrics

Behaviour of a plain knit fabric

The behaviour of a knitted fabric can be divided into two main steps separated by a jamming transition phase (Fig. 4a). The first step involves a slipping of the yarns at the contact point between adjacent loops. A change of the loop shape occurs by yarn bending and results in a linear behaviour with an important deformation of the fabric under low loads. This phenomenon occurs until the jamming transition phase, in which a transverse compression of the yarns occurs due to fibres rearrangement at the cross-over. The yarns bending tends toward a limit while the tensile load increases. At the end of this phase, the stitch is fully jammed. The load is then carried out by the stitch legs, subjected to tensile stresses, whereas shear stresses increase in the cross-over areas. A fast increase of the load is thus observed, resulting in a second almost linear domain of the load-strain curve. The first step is characterized by a critical strain (e c ) which corresponds to the strain at the intersection of both linear domains. The last step is characterized by the fabric stiffness after full jamming of the stitches, corresponding to the slope of the curve, and by the failure load.

Fig. 4b illustrates the tensile behaviour of the plain-knit according to the loading direction and the fibre type fabrics (E-glass, basalt and carbon). Table 4 reports the values of the characteristic parameters of the tensile behaviour. Irrespective of the fibre type and loading direction, plain-knit fabrics behave similarly. A higher deformation ability and lower mechanical properties are observed in the course-direction. These properties result from the loop geometry related to the knitting process: the yarns run along the courses and are interconnected in the wale-direction. In Fig. 5 are shown schematic views of the fabric structure at initial state and stretched in wale or course-directions. Stretched in the course-direction, the stitches lose their loop shape to take a curvy form, whereas in the wale-direction, the stitches present a stretched loop shape at jamming. In the first step, considering the low loads involved, the fibre type has a negligible effect on the fabric macroscopic behaviour and the critical strains reached are similar for a given loading direction (Table 4). After jamming, one yarn per course is subjected to tensile stress while two yarns per wale are loaded. It explains the higher rigidity and the higher strength reached by plain-knits loaded in the wale-direction. In this part, the fibre type is the most influent parameter. Best properties are obtained with carbon fibres whereas glass and basalt plain-knits exhibit similar properties despite the slightly higher values reached with basalt fibres.

Inlay yarn fabric behaviour

The tensile behaviour of plain-knits with inlay yarns is plotted in Fig. 6 and the properties are summarized in Table 4. In the wale-direction, the insertion of inlays tends to decrease the strain at jamming, but this effect is slight compared to the standard deviation of the experimental results. Each yarn diameter being similar to each others, the strain levels reached are of the same order of magnitude, irrespective of the fibre type. Fig. 5 depicts a schematic view of a plain-knit with inlay yarns stretched in the wale-direction. In this case the inlay yarns are not tightened and do not contribute to the fabric behaviour except by slightly modifying the size of the stitches. P1, P2 and P3 fabrics behave similarly to P0. The differences observed are attributed to variations in the stitches dimensions and to inlay yarns that decrease slightly the critical strain. Considering the experimental dispersion, it could be concluded that the mechanical properties are almost the same with glass and basalt fibres. With carbon fibres, the amount of broken fibres increases with the number of inlay yarns, resulting in a loss of properties.

In the course-direction, the pneumatic grips maximal load limit was reached before P2 and P3 knits failure. The inlay yarns control the behaviour of the fabric and the plain-knit performs the function of carrying the inlays. Thus, the critical strain of the fabric is the strain required to tighten the floated yarns (Fig. 6), decreasing sharply the critical strain e c (Fig. 6). The fibre orientation of the straightened yarns close to the loading direction has for consequence a high increase of the fabric rigidity. The curves related to P1 fabrics show a second increase of the load after the rupture of the inlay yarns: plain stitches are loaded and the fabric behaves like a plain-knit stressed in the course-direction, but at lower strains levels due to the inlay yarns. The fabric rigidity increases systematically with the number of inlays, however this increase is more or less important depending on the fibre type. Comparing P0 and P3, an improvement of the rigidity of 770% is obtained with glass fibres whereas it is only of 308% with carbon fibres and of 502% with basalt fibres. Although carbon fibres are the most efficient, the properties obtained with the glass fibres are the best. Similarly, between P0 and P1, the ultimate strength has been multiplied by 5 with glass and basalt fibres and only by 2.5 with carbon fibres. Indeed, it is not possible to ensure a perfectly identical length of the different inlay yarns, preventing them to be tighten simultaneously. These trends are also attributed to the amount of broken fibres in the fabric: producing more complex fabrics induces an increase of damaged fibres. The shear strength depends on the filament cross-section and seems to be the key parameter that controls the ability of the yarn to be knitted without damage.

Behaviour of knit fabric reinforced composites

Specimen and testing method

Tensile tests were performed in the directions 0°(wale-direction), 45°and 90°(course-direction) on 20 Â 180 mm 2 plain-knit (P0) and on 25 Â 200 mm 2 P1, P2, and P3 composite specimens according to ASTM D3039/D3039M-00 recommendations. P1, P2 and P3 composite specimens are larger in order to get a full representative elementary cell of the knitted reinforcement (see Fig. 2.) in the specimen width. Four 20 Â 40 or 25 Â 40 mm 2 aluminium tabs were bonded on each specimen in order to limit stress concentration in the clamping areas. The tensile tests were carried out on an INSTRON 5800R electromechanical testing machine at the constant crosshead displacement rate of 0.5 mm min À1 . The longitudinal strain is measured with a 50 mm gauge length extensometer. Each configuration was tested three times. In this section a representative test among the three was selected to plot and discuss the results.

Mechanical testing results

Due to the matrix involved in this study, characterised by a low strain at break, the fibrous architecture does not change during the tests, and the composite behaviour is mainly related to the initial shape of the knitted loop.

The stress-strain curves are plotted in Fig. 7 and the mean values of the mechanical properties are given in Table 5. The behaviour of knit reinforced composites is elastic damageable with various behaviours in the non-linear domain. They exhibit a high behaviour diversity depending on the loading direction, fibrous This figure was drawn using pictures adapted from [START_REF] Araújo | Modelling and simulation of the mechanical behaviour of weft knitted fabrics for technical applications -Part 2: 3D model based on the elastica theory[END_REF][START_REF] Bekisli | Mechanical properties of thermoformed structures with knitted reinforcement[END_REF].

architecture and fibre type. In all cases the fibre reinforcement increases the material elastic modulus when comparing to the net resin sample behaviour. P0 reinforced composites exhibit an anisotropic behaviour with higher properties in 0°direction and similar properties in 45°and 90°directions. In this latter direction, the ultimate strength is lower than the net resin strength. The behaviour of the P0 composites reinforced with basalt and carbon is close to a brittle elastic one, while the P0 composite reinforced with glass fibres, whose properties are the lowest, exhibit a wider non-linear domain before failure (45°and 90°).

Inserting one inlay yarn per course in plain-knit reinforcement has for consequence a slight decrease of the composite properties in the 0°direction attributed to a different arrangement of the fibres distribution in the material at similar fibre volume fractions. The insertion of one inlay yarn has no significant effect in the directions 45°and 90°. A slight improvement is obtained with glass fibres whereas a small decrease of modulus and strength is observed in the direction 45°with basalt and carbon fibres. Similar properties are obtained in the direction 90°except a wider non-linear domain and a slight increase of strength with carbon fibres. Indeed, the floated yarns are not straight but curved, limiting a major orientation of the fibres in the direction 90°, and limiting the expected improvement of the properties in this direction. In addition, a higher porosity is observed in P1 (Table 3).

Increasing the number of inlay yarns per course leads to an important improvement of the properties in the direction 45°a nd an even larger increase in the direction 90°whereas no significant change is observed in the direction 0°. However, in the loading direction 0°, after a loss of properties with the insertion of the first inlay yarn (P1), an increase of the properties is observed with the insertion of the second inlay yarn (P2), attributed to the increase of the fibre content, and again a decrease in properties is observed with the third inlay yarn (P3). In the others loading directions, between P0 and P3, according to the fibre type, an increase of 60-105% of the modulus and 61-134% of the ultimate stress is experienced in the testing direction 45°, and an increase of 76-134% of the modulus and 212-288% of the ultimate stress is reached in the loading direction 90°. This is due to a cross-effect of fibre content and fibre distribution. In a general way, the stressstrain curves show that the non-linear domain of composites reinforced by plain-knits with inlay yarns is relatively strong in comparison to P0 reinforced composites. Increasing the number of floated yarns leads to more diffuse damages and to higher strain at break. The effect of the fibre type depends on the fibrous architecture and the loading direction. In the direction 0°the fibre type effect on the properties decrease is similar irrespective of the fibre types: it is significant for the first inlay yarn whereas a low decrease is noticed when increasing the number of inlay yarns. On the contrary, the properties in the directions 45°and 90°are directly related to the number of inlays. In these directions, the more the number of inlay yarns is, the higher the properties of carbon fibres increase, as shown by the increasing gap between the results obtained with carbon fibres and those obtained with glass and basalt fibres. In all cases, best properties are obtained with carbon fibres and this is in agreement with the fibre properties given in Table 1.

Fracture analysis

Figs. 8 and9 show photographs and images of cracked specimens as well as a schematic view of the fracture surface on the knit samples. In the loading directions 0°and 90°(Fig. 8a andb), the failure occurs systematically in a plane perpendicular to the loading direction. For P0 composite specimens tested in the direction 0°, the cracks propagate along a crenel-shape pattern in the vicinity of the fracture (Fig. 8a; ESEM image 1). The crack initiates perpendicularly to the loading direction by a separation between a fibre tow and the matrix and propagates between a fibre-rich area and a resin-rich area up to the cross-over points. When the crack meet the cross-over area, a bifurcation occurs and the material damages under shear stresses up to the leg of the stitch. The energy required by this mechanism allow the same phenomenon to take place elsewhere before to go on. At least one crack propagates along adjacent stitches until the specimen failure throughout the width of the sample. The two parts of the specimen are then connected by strands of fibres that are then loaded up to the failure (Fig. 8a; SEM image 2). This explains the high strain values reached with P0 composites in the direction 0°in contrast to 45°and 90°( Fig. 7).

When inserting inlay yarns in the reinforcement, the rupture plane remains transverse to the loading direction between the rows of stitches (Fig. 8c). Whatever the fibre type or the number of inlays, the composite damages similarly. The cracks are only located in the failure area. The cracks initiate between resin-rich and fibre-rich areas and propagate along floated yarns. As for P0 composites, the two parts of the specimen are then connected by strands of fibres of the plain knit stitches which are then loaded up to the failure.

In the loading direction 90°, the fracture occurs between resinrich and fibre-rich areas, by debonding between fibre tows and matrix whatever the fibrous architecture (Fig. 8b andd). The cracks follow perfectly the stitches legs of a single wale. The inlay yarns plays the role of slowing down the cracks propagation by increasing the fibre volume fraction perpendicular to the load in this area, which leads to increase the material strength.

In the P0 specimens tested in the direction 45°, the fracture occurs in the same way as for those tested in the direction 90°(between a rich and a poor fibre region). A 45°oriented fracture is then obtained (Fig. 9a), and this explains the equivalent properties obtained with P0 composites tested at 45°and 90°. When inserting one inlay yarn, the behaviour is quite similar (Fig. 9b). By inserting two inlay yarns, two cases may arise (Fig. 9c). In the first case, the crack propagation is similar to those of the composites tested in the direction 90°(perpendicular to the inlay yarns, Fig. 8d). This is observed with carbon fibres reinforced composites. In the second case, observed with basalt and glass fibres composites, the crack propagation is similar to the P2 composites tested in the direction 0°(Fig. 8c), i.e. along the floated yarns, leading to a 45°oriented fracture surface occurring along the inlay yarns. This latter failure mode systematically occurs with P3 specimens tested in the direction 45°. With one inlay yarn, the cross-over area is strongly sheared, especially in the tuck stitch area since the floated yarn strengthens the material. With three inlays, the floated yarns constitute a continuous unidirectional stripe. High shear strains arise thus in the vicinity of these stripes in areas where the fibres are mainly wale-wise oriented. The failure occurs between these two areas. With two inlays, the tucked wale is not efficiently reinforced by the floated yarns since they are not perfectly parallel to the course (see Fig. 2). Both failure mechanisms are thus possible depending on the fibres properties and on the amount of broken fibres.

Conclusion

Insertion of up to three inlays yarns in a plain knit fibre architecture was investigated in order to improve the mechanical behaviour of knit-reinforced composites. Each knit fabric was manufactured with E-glass, basalt and carbon fibres, using the same stitch size and yarn diameter and with a tuck every four wale for the inlaid knits.

The tensile behaviour of each knit in the wale-direction shows that the wale-wise behaviour is not significantly affected by the inlays, except by modifying the cross-over areas with the tucks. At the contrary, the course-wise direction is significantly reinforced and about five times less deformable with inlay yarns. Indeed, the behaviour of the knit in this case is mainly the behaviour of the inlays. The effect of the fibre type on the behaviour of the knit is negligible compared to the effect of the inlays.

The stress-strain behaviour of the knit reinforced composite material depends on the fibre type but also strongly on the knit construction. Carbon fibres provide the best mechanical properties but the results are not proportional to the intrinsic properties of the fibres when comparing with E-glass and basalt. Inlaying floated yarns does not allow to improve the mechanical behaviour of the composite since it decreases the strength in the direction 0°, but it makes possible to control the anisotropy of the material properties. Highest mechanical properties are reached in the direction 0°f or plain-knit reinforced composites, in 90°direction for P3 knitted composites and an almost isotropic material is obtained for the P2 knitted composites.

Fig. 1 .

 1 Fig. 1. Plain-knit basic architecture (P0): (a) photograph of a basalt fabric; (b) testing directions and (c) representative elementary cell (REC) and geometrical characteristics of the loop.

3. 1 .Fig. 2 .

 12 Fig. 2. Photographs, schematic views and representative elementary cell (REC) of the knitted fabrics: (a) with 1 inlay yarn (P1); (b) with 2 inlay yarns (P2) and (c) with 3 inlay yarns (P3).

Fig. 4 .

 4 Fig. 4. General behaviour of a plain knit fabric: (a) characteristic parameters, and (b) tensile load-true strain curves of E-glass (G), basalt (B) and carbon (C) plain-knit fabrics in wale and course-directions.

Fig. 5 .

 5 Fig. 5. Scheme of a plain-knit fabric and plain-knit with inlay yarns at (a and d) initial state, (b and e) stretched in wale direction and (c and f) stretched in course direction.This figure was drawn using pictures adapted from[START_REF] Araújo | Modelling and simulation of the mechanical behaviour of weft knitted fabrics for technical applications -Part 2: 3D model based on the elastica theory[END_REF][START_REF] Bekisli | Mechanical properties of thermoformed structures with knitted reinforcement[END_REF].

Fig. 6 .ig. 7 .

 67 Fig. 6. Wale and course direction load versus true strain curves for dry glass, basalt and carbon fibre knitted fabrics. The dotted lines for P2 and P3 knits in the course direction indicate that the failure has not been reached.

Fig. 8 .

 8 Fig. 8. Failure of plain-knit with and without inlay yarns reinforced composites stressed in 0°(wale) and 90°(course) directions.

Table 1

 1 Fibre, yarn and resin properties.

	Fibre type	E-glass	Basalt	Carbon	Resin type
		Owens corning Vetrotex T30 111AX23 Basaltex BCF Tohotenax HTA 5131 6 K Z10 Thermoset epoxy resin SICOMIN SR1710 + SD8822
					hardener
	Fibre properties				Tensile properties
	Filament diameter	16 lm	1 3 lm	7 lm	Young's modulus	3680 MPa
	Fibre density	2620 kg m À3	2700 kg m À3 1760 kg m À3	Rupture stress	85 MPa
	Young's modulus	80 GPa	84 GPa	238 GPa	Strain at break	3.1%
	Yarn properties				Bending properties
	Number of filaments 1200	1672	6000	Modulus	3720 MPa
	Yarn diameter	0.58 mm	0.56 mm	0.57 mm	Maximal stress	136 MPa
	Linear density (Dy)	600 tex	600 tex	400 tex	Strain at break	5.20%
	Bending stiffness	0.31 N mm 2	0.20 N mm 2	0.16 N mm 2

Table 2

 2 Knitted fabrics characteristics.

	P0	P1	P2	P3
	Thx (mm)

E-glass 1.60 ± 0.00 1.92 ± 0.06 2.22 ± 0.13 2.82 ± 0.75 Fabric thickness (measured under 0.005 MPa pressure) Basalt 1.55 ± 0.09 1.83 ± 0.06 2.20 ± 0.05 2.78 ± 0.75 Carbon 1.53 ± 0.03 1.80 ± 0.05 2.08 ± 0.03 2.62 ± 0.60 M S (g m À2 ) E-glass 792 ± 60 1028 ± 181 1382 ± 167 1589 ± 160 Areal density Basalt 797 ± 113 1044 ± 227 1337 ± 161 1572 ± 173 Carbon 566 ± 32 664 ± 136 874 ± 152 1072 ± 108 Upper mould Preform Vacuum outlets Injection intlet Resin tank Lower mould Spacer Fig. 3. VARTM processing conditions.

Table 3

 3 Composite plate characteristics.

		Reinforcement	P0	P1	P2	P3
		Number of plies	3	2	2	2
	Plate thickness (mm)	E-glass	3.15	3.18	3.19	3.22
		Basalt	3.15	3.17	3.19	3.21
		Carbon	3.15	3.18	3.17	3.17
	Fibre content (%)	E-glass	31	29	38	43
		Basalt	32	31	35	41
		Carbon	32	30	37	42
	Porosity (%)	E-glass	2.8	5.4	7.1	5.4
		Basalt	3.4	8.4	7.5	6.6
		Carbon	2.6	10.1	4.1	3.2

Table 4

 4 Tensile properties of the knitted fabrics.

	Knitted fabric	Fibre type	Critical stretched strain e c	Fabric rigidity (N/mm/mm)	Ultimate tensile load Fr (N)
			Wale	Course	Wale	Course	Wale	Course
	P0	Glass	0.57 ± 0.02	0.81 ± 0.03	1450 ± 63	443 ± 27	57.2 ± 7.5	18.3 ± 1.9
		Basalt	0.59 ± 0.02	0.91 ± 0.01	1573 ± 384	482 ± 71	70.0 ± 8.0	21.5 ± 3.1
		Carbon	0.52 ± 0.03	0.88 ± 0.02	2292 ± 351	768 ± 93	111.0 ± 3.0	37.2 ± 3.9
	P1	Glass	0.58 ± 0.02	0.19 ± 0.01	1740 ± 138	1917 ± 176	54.0 ± 1.6	90.9 ± 12.0
		Basalt	0.54 ± 0.01	0.15 ± 0.01	1866 ± 164	2315 ± 94	57.8 ± 1.9	108.6 ± 5.5
		Carbon	0.50 ± 0.01	0.15 ± 0.05	1560 ± 52	2349 ± 366	56.2 ± 1.8	97.4 ± 8.6
	P2	Glass	0.54 ± 0.01	0.18 ± 0.01	1525 ± 150	3513 ± 528	52.2 ± 2.6	-a
		Basalt	0.52 ± 0.02	0.16 ± 0.02	1542 ± 57	2744 ± 177	58.1 ± 3.6	-a
		Carbon	0.51 ± 0.01	0.16 ± 0.02	1536 ± 74	2422 ± 142	64.1 ± 0.7	-a
	P3	Glass	0.42 ± 0.03	0.13 ± 0.01	1538 ± 121	3856 ± 269	66.2 ± 3.9	-a
		Basalt	0.43 ± 0.03	0.09 ± 0.01	1756 ± 165	2902 ± 86	72.0 ± 3.3	-a
		Carbon	0.47 ± 0.00	0.13 ± 0.01	1370 ± 217	3137 ± 105	71.5 ± 3.8	-a

a Allowable load of the pneumatic grips reached before fabric break.

Table 5

 5 Tensile properties of the knit reinforced composites.

	Fibre type	Knitted fabric	Longitudinal elastic modulus (GPa)		Ultimate stress (MPa)
			0°45°90°0°45°90°G
	lass	P0	11.2 ± 0.1	8.1 ± 0.2	8.2 ± 0.0	112.2 ± 4.8	54.1 ± 2.2	42.3 ± 3.4
		P1	9.8 ± 0.3	9.0 ± 0.4	8.8 ± 0.2	90.7 ± 3.2	66.4 ± 0.2	57.1 ± 3.0
		P2	10.8 ± 0.2	11.3 ± 0.3	12.2 ± 0.2	102.8 ± 2.4	86.7 ± 2.1	124.1 ± 5.1
		P3	12.1 ± 0.2	13.3 ± 0.4	17.2 ± 0.6	82.4 ± 0.5	87.3 ± 2.7	164.2 ± 1, 2
	Basalt	P0	12.8 ± 0.1	8.5 ± 0.1	9.3 ± 0.8	136.9 ± 5.4	54.8 ± 0.9	51.1 ± 0.3
		P1	9.7 ± 0.4	8.3 ± 0.3	9.0 ± 0.2	84.2 ± 9.2	47.7 ± 1.0	52.5 ± 2.6
		P2	10.8 ± 0.2	10.8 ± 0.2	12.0 ± 0.8	97.1 ± 3.7	85.8 ± 4.5	130.0 ± 6.6
		P3	11.8 ± 0.3	13.5 ± 0.1	16.4 ± 0.8	87.9 ± 4.6	95.8 ± 6.0	160.9 ± 5.4
	Carbon	P0	17.5 ± 0.7	8.9 ± 0.0	9.3 ± 0.7	152.1 ± 0.2	55.6 ± 3.2	63.3 ± 3.6
		P1	12.7 ± 0.1	8.8 ± 0.2	10.6 ± 1.3	107.9 ± 15	46.0 ± 0.7	75.2 ± 7.6
		P2	13.7 ± 0.7	15.0 ± 0.8	19.2 ± 1.4	121.5 ± 4.6	106.6 ± 6.8	163.4 ± 6.3
		P3	12.7 ± 0.8	18.2 ± 1.0	21.8 ± 0.2	104.9 ± 3.0	130.2 ± 5.2	197.4 ± 0.0
	Resin		3.62			65.5

Acknowledgements

This study was carried out with the financial support of the Region Midi-Pyrénées in the framework of the project CORTEX. The authors would like to thank JTT Composite, and especially MM. Nicolas Dumont and Francis Planel for the design and manufacturing of the knits. wale wale wale wale Fig. 9. Failure of plain-knit with and without inlay yarns reinforced composites stressed in 45°direction.


Antoine Pierquin

Stephane Brisset

T Henneron

Multidisciplinary optimization formulation for the optimization of multirate systems

Keywords: Electromagnetic coupling, multi-physics and coupled problems, waveform relaxation method, multidisciplinary optimization

come

Multidisciplinary optimization formulation for the optimization of multirate systems

I. INTRODUCTION

M ULTIDISCIPLINARY optimization (MDO) considers the different ways to insert and use a multiphysic model into an optimization problem. Indeed, the model can be more or less integrated to the optimization process by changing the formulation of the optimization problem. Different levels of integration are possible when the multiphysic modeling is based on a fixed-point approach. In this case, each subsystem (physical phenomenon, device or discipline) of the multiphysic model is solved with an adapted numerical model. Then, the coupling is done by using a solution from a subsystem as a source for another subsystem. Lastly, the consistency of the coupling between the subsystems is guaranteed in an iterative way. In optimization, the fixed-point approach can be exploited by adding the fixed-point as an integral part of the optimization problem.

These considerations are known and used since a long time in optimization in the static case [START_REF] Braun | Development and application of the collaborative optimization architecture in a multidisciplinary design environment[END_REF]- [START_REF] Ayed | Multidisciplinary optimization formulations benefits on space mapping techniques[END_REF]. Two main approaches are available: the multidisciplinary feasibility (MDF) and the individual disciplinary feasibility (IDF) [START_REF] Cramer | Problem formulation for multidisciplinary optimization[END_REF]. In the first one, the fixed-point loop is performed systematically when the multiphysic model is evaluated. In the second one, the fixed-point criterion is added as a constraint of the problem. Consequently, with the MDF approach, the consistency of the coupling is guaranteed at each evaluation of the model whereas with the IDF approach, the consistency of the model is only guaranteed for the optimal solution while the number of model evaluations is reduced.

In the dynamic case, fixed-point approach also exists, known as waveform relaxation method (WRM) [START_REF] Lelarasmee | The waveform relaxation method for time-domain analysis of large scale integrated circuits[END_REF], [START_REF] Bartel | Dynamic iteration for coupled problems of electric circuits and distributed devices[END_REF] that is welladapted for multirate system, a system whose subsystems have different time-constants. However, to the best of our knowledge, MDO strategies have not been studied in dynamic case. Authors already studied the MDF [START_REF] Pierquin | Benefits of waveform relaxation method and output space mapping for the optimization of multirate systems[END_REF] and IDF approaches but the use of the last one was problematic for gradient-based optimization methods. Indeed, the computation of the gradient was too time-consuming because it was computed by finite differences and the number of optimisation variables is high with IDF. This paper demonstrates that the partial derivatives of the objective and constraints versus all time-dependant variables can be obtained by using only one additional evaluation of the model. Thus, the computation of the gradient requires a small number of model evaluations, allowing a reduction of the optimization time compared to the MDF approach.

In a first time, the WRM will be presented, then the MDF and IDF strategies in dynamic case. Lastly, they will be applied on an academic example and analysed in terms of complexity and computation time.

II. WAVEFORM RELAXATION METHOD

The global system E ( ẇ(t), w(t)) = 0 is considered to be decomposed in r subsystems such that each subsytem verifies

E j ( ẇj (t), w j (t), d j (t)) = 0, (1)
with d j (t) = [w 1 (t), . . . , w j-1 (t), w j+1 (t), . . . , w r (t)].

The WRM, also called dynamic iteration, solves iteratively the r subsystems of differential algebraic equations to produce an approximation of the exact solution. At the k th iteration and for the subsystem j, the following problem is solved:

E k j ẇk j (t), w k j (t), d k j (t) = 0, (2)
with d k j (t) = w k 1 (t), . . . , w k j-1 (t), w k-1 j+1 (t), . . . , w k-1 r (t) with a Gauss-Seidel scheme (see [START_REF] Lelarasmee | The waveform relaxation method for time-domain analysis of large scale integrated circuits[END_REF] for other schemes). The solution w k is obtained from the previous solution w k-1 , so the solution at the k th iteration can be written by introducing an operator Ψ as and for each subsystem w k j = Ψ j (w k-1 ), j = 1 to r. Note that if w k-1 = w k , then w k verifies (1) and is solution of the global system. So the opterator Ψ is a fixed-point operator [START_REF] Bartel | Dynamic iteration for coupled problems of electric circuits and distributed devices[END_REF]. Finally, the algorithm stops when w k (t) is close enough to w k-1 (t). At the end of the iterative process, the solution is w K (t) = Ψ K w 0 (t) .

w k = Ψ(w k-1 ), (3)
To apply a gradient-based algorithm to solve the fixed-point problem, or later in optimization in the IDF formulation, the Ψ-operator will be approximated in the neighbourhood of w(t) by the first order Taylor expansion

Ψ(w(t)) Ψ( w(t)) + ∇Ψ.(w(t) -w(t)). (4)

III. MULTIDISCIPLINARY OPTIMIZATION A. Multidisciplinary feasibility

Let us consider the optimization problem

x = arg min x f (x) such that g(x) ≤ 0, (5)
and

h(x) = 0, (6)
with x the degrees of freedom, f the objective function and g and h the constraints. We suppose that to obtain the output g(x) and h(x), a system modeled by WRM is solved. With an MDF approach, at each evaluation of g, the iterative procedure of WRM is performed until the fixed-point criterion is verified. Thus, the consistency of the coupling is insured for all evaluations of the model. The Fig. 1 shows the principle with two subsystems Σ 1 and Σ 2 : the waveform solution w 1 of Σ 1 is the source of Σ 2 , and reciprocally the waveform solution w 2 of Σ 2 is the source of Σ 1 .

If the operator Ψ defined in (3) is applied K times on average, the number of evaluations of Ψ is approximately n eval ×K at the end of the optimization process, with n eval the number of model evaluations during the optimization process.

Subsystem

Subsystem

w 1 =Ψ 2 (w 2 ) w 2 =Ψ 1 (w 1 ) Optimiser x f (x ), g (x ), h (x ) Ψ K (w ) Σ 1 Σ 2
Fig. 1. Scheme of the MDF principle.

B. Individual discipline feasibility

With an IDF approach, the initial problem ( 6) is modified. The discretized waveform w(t) is added to the optimization variables and the fixed-point condition is also added as a constraint of the problem. Only one evaluation of Ψ is done per model resolution, the consistency of the coupling being guaranteed at the end of the optimization process. The optimization problem to solve is

[x , w ] =arg min x,w f (x, w) such that (7) g(x, w) ≤ 0, (8) h(x, w) = 0, (9)
w T -Ψ (w) = 0. ( 10
)
One negative point of the IDF approach is to considerably increase the number of optimization variables. The size of the problem can become too high to be solved by non-linear optimization methods.

The Fig. 2 shows the principle with the same two subsystems Σ 1 and Σ 2 of Fig. 1, in two possible configurations. The subsystems are solved in parallel from the input given by the optimiser in Fig. 2(a); the subsystems are solved sequentially in Fig. 2(b). With this second configuration, the waveform w added as optimization variables can be reduced to w 1 .

At the end of the optimization process, the number of evaluations of Ψ is ñeval , with ñeval the number of model evaluations. The optimization duration will be reduced if ñeval < n eval × K.

Optimiser Subsystem Subsystem Ψ 1 (w 1 ) Ψ 2 (w 2 ) w 1 w 2 x , w f (x ) g (x ) h (x ) Ψ (w )-w Σ 1 Σ 2
(a) IDF principle with parallel solving of the subsystems.

Optimiser Subsystem Subsystem Ψ 1 (w 1 ) Ψ 2 (w 2 ) w 1 w 2 x , w f (x ) g (x ) h (x ) Ψ (w )-w Σ 1

IV. APPLICATION EXAMPLE

A. Problem formulations

A LC filter supplying a transformer is considered. The voltage source v s is a 50 Hz fundamental signal and a 20 kHz sine perturbation. Different time-steps can be used to model this device by WRM (Fig. 3), with a 2D finite element model (FEM) for the simulation of the transformer. The time-step into the FEM is fitted to the 50 Hz dynamic, whereas in the circuit part, it is fitted to the 20 kHz dynamic. With a Gauss-Seidel scheme, the circuit model and the FEM are solved sequentially. At each WRM iteration, a current waveform i k-1 (t) is imposed as a source into the circuit model. Its resolution gives the voltage waveform v k (t) imposed to the transformer, then the resolution of the FEM gives the current i k (t) that will be the current source of the circuit model at the next iteration. The time-discretisation being different, linear interpolation is used to obtain the waveform for all time t.

At each iteration i k = Ψ(i k-1 ), with i k = [i k (t 1 ), i k (t 2 ), . . . , i k (t N )] T and t n , n = 1 to N the time discretisation, and so the algorithm produces i K = Ψ K (i 0 ). 3. Device split to apply the WRM.

v i i L i v v s i v Fig.
The optimization aims to minimize the mass of the transformer by acting on the width L and the height H of the transformer (Fig. 4). Moreover, the root mean square current i rms into the transformer has to be equal to 3 A. The initial problem, solved with an MDF approach, is

             min H,L m(H, L), 20cm ≤ H ≤ 40cm, 12cm ≤ L ≤ 24cm, H -2L 3 > 0, i rms = 3 A. (11)
The sequential solving of the models (circuit then FE) allows an IDF formulation with only one waveform added as optimization variables, like in Fig. 2(b). The discretized current waveform i is added to the optimization variables. The sequential solving of the circuit and of the FEM produces the output current waveform i out = Ψ(i): the fixed-point criterion i out = i is added as a constraint into the optimization problem. The waveform of the voltage v is juste an intermediary variable that does not take part of the optimization problem. By adding only i to the problem, it avoids increasing too much the dimension of the optimization problem. We note i n = i(t n ) and i out n = i out (t n ). With these notations, the problem to solve with the IDF formulation is

                 min H,L,i m(H, L), 20cm ≤ H ≤ 40cm, 12cm ≤ L ≤ 24cm, H -2L 3 > 0, i rms = 3 A, i n -i out n = 0, ∀n = 1, . . . , N. (12)

B. Computation of the jacobian

The discretized waveform of the current is added to the optimization variables and the constraints i n -i out n = 0, for n = 1 to N ensure i -Ψ(i) = 0, namely the consistency of the coupling. One difficulty is to obtain the jacobian of the operator Ψ to use a gradient-based optimization algorithm. Most of the time, the gradient is computed by finite differences for each optimization variable. Nevertheless, this is not required for the gradient with respect to i n Due to the implicit Euler scheme used to solve the differential equations, the jacobian ∇Ψ is a triangular matrix if the time-discretisation is the same for all the subsystems. Indeed, ∂i out q ∂ip = 0 if q < p, because a perturbation at the moment t p in input does not have any influence on the output for the previous moments. Moreover, we have ∂i out p ∂ip = ∂i out q ∂iq , p = q and consequently ∂i out p+ ∂ip = ∂i out q+ ∂iq , ∀ ≥ 0. This means that the matrix ∇Ψ has identical values onto its diagonals, and only ∂i out ∂i1 is needed to complete the matrix. . . .

J =         
∂i out 1 ∂i1          (13)
To explain the last point, let us consider more generally the WRM applied with a Gauss-Seidel scheme to the following linear differential systems

∂ ∂t X 1 = A 1 X 1 + B 1 X in (14) ∂ ∂t X out = A 2 X 1 + B 2 X out . ( 15
)
where A 1 , A 2 , B 1 , B 2 are constant matrix. With an Euler scheme in time, the systems become

X 1 n = 1 δt I 1 -A 1 -1 1 δt X 1 n-1 + B 1 X in n ( 16
)
X out n = 1 δt I 2 -B 2 -1 1 δt X out n-1 + A 2 X 1 n , (17)
with I 1 , I 2 the identities matrix with the adadted size.

The derivative of the previous expressions with respect to

X in n gives ∂X 1 n ∂X in n = 1 δt I 1 -A 1 -1 B 1 (18) ∂X out n ∂X in n = 1 δt I 2 -B 2 -1 A 2 ∂X 1 n ∂X in n , (19)
and so

∂X out n ∂X in n = 1 δt I -B 2 -1 A 2 1 δt I 1 -A 1 -1 B 1 , ∀n.
(20) This proves that ∂i out p ∂ip = ∂i out q ∂iq , p = q and consequently ∂i out p+ ∂ip = ∂i out q+ ∂iq , ∀ ≥ 0 by induction. This can be proved with more than 2 subsystems if the solving of the subsystems is totally sequential, when the solution of a subsystem is used as a source for the next one.

Finally, the computation of

∂i out
∂i1 is obtained for all by only one calculation by finite differences, and allows to have the complete matrix ∇Ψ. The application Ψ is evaluated with the waveform input i:

i out = Ψ(i). (21)
Then we subject the input i 1 to a small perturbation δi, before to evaluate Ψ with the perturbed waveform:

i out δ = Ψ(i + δie 1 ), (22)
with e 1 the canonical basis. Lastly, the gradient can be approximated ∂i ∂i

1 i out δ -i out δi . (23)

C. Results

Optimizations using the sequential quadratic programming (SQP) method of Matlab R software are done for the MDF and IDF formulations. SQP finds local optimum and needs the computation of the gradient of the objective and constraint functions. The gradient of the constraints is computed by finite differences with the technique presented in section IV-B for ∇Ψ. So at each call of the constraint function, the application Ψ is evaluated four times: one times to have the value of the function at a given point, then three times to obtain the derivative with respect to H, L and ∇Ψ.

To compare the two approaches, 10 random initial points are generated. For 8 points, the IDF formulation converges to the same solution than the MDF one. For these 8 points, the average relative error is 0.0012%. But with the IDF approach, the number of evaluations of the Ψ operator is reduced (Fig. 5). Each evaluation of Ψ corresponds to the solving of a FE problem in the time-domain, that is time-consuming. So the reduction of the number of FE problem solvings leads most of the time to a decrease of the optimization time. Nevertheless, the optimization process can evaluate a lot of times the objective function without computing the constraints. This explains that the optimization duration stays superior with the IDF formulation for the 7 th point. Finally, the optimization duration is globally lower with the IDF approach. The speedup factor is 3 on average, and even greater than 7 for one point (Fig. 6).

V. CONCLUSION

The WRM used to co-simulate multirate systems has been added into the optimization problem thanks to the IDF formulation. The dimension of the problem is increased but the number of evaluations of the model decreases through an efficient way to compute the derivatives of time-dependent variables. This leads to a significant reduction of the optimization time.

Σ 2 (Fig. 2 .

 22 Fig. 2. Scheme of the IDF principle.

3 Fig. 4 .

 34 Fig. 4. Geometry of the transformer.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Number of evaluations of Ψ operator.

  0018-9464 c 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications standards/publications/rights/index.html for more information. (Inserted by IEEE.)

Manuscript received July 5, 2015; revised May 15, 2015 and June 1, 2015; accepted July 1, 2015. Date of publication July 10, 2015; date of current version July 31, 2015. Corresponding author: A. Pierquin (e-mail: ant.pierquin@gmail.com) Digital Object Identifier (inserted by IEEE).