<figure><img src="image_1.png" /></figure>

HAL
open science

# On the use of PGD for optimal control applied to automated fibre placement

Nicolas Bur, Pierre Joyot

► To cite this version:

Nicolas Bur, Pierre Joyot. On the use of PGD for optimal control applied to automated fibre placement. 20th International Esaform Conference on Material Forming: Esaform 2017, Apr 2017, Dublin, Ireland. pp.170006, 10.1063/1.5008204. hal-01736984

HAL Id: hal-01736984

https://hal.science/hal-01736984

Submitted on 22 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the Use of PGD for Optimal Control Applied to
Automated Fibre Placement

N. Bur<sup>1,a)</sup> and P. Joyot<sup>1,2,b)</sup>

<sup>1</sup>ESTIA-Recherche, Technopôle Izarbel, 64210 Bidart, France

<sup>2</sup>I2M, UMR CNRS 5295, 33405 Talence Cedex, France

**Abstract.** Automated Fibre Placement (AFP) is an incipient manufacturing process for composite structures. Despite its conceptual simplicity it involves many complexities related to the necessity of melting the thermoplastic at the interface tape-substrate, ensuring the consolidation that needs the diffusion of molecules and control the residual stresses installation responsible of the residual deformations of the formed parts.

The optimisation of the process and the determination of the process window cannot be achieved in a traditional way since it requires a plethora of trials/errors or numerical simulations, because there are many parameters involved in the characterisation of the material and the process.

Using reduced order modelling such as the so called *Proper Generalised Decomposition* method, allows the construction of multi-parametric solution taking into account many parameters. This leads to *virtual charts* that can be explored on-line in real time in order to perform process optimisation or on-line simulation-based control. Thus, for a given set of parameters, determining the power leading to an optimal temperature becomes easy.

However, instead of controlling the power knowing the temperature field by particularizing an abacus, we propose here an approach based on optimal control: we solve by PGD a dual problem from heat equation and optimality criteria. To circumvent numerical issue due to ill-conditioned system, we propose an algorithm based on Uzawa's method. That way, we are able to solve the dual problem, setting the desired state as an extra-coordinate in the PGD framework. In a single computation, we get both the temperature field and the required heat flux to reach a parametric optimal temperature on a given zone.

## INTRODUCTION

Automated Fibre Placement (AFP) is one the main technologies employed today to manufacture advanced composite laminates from unidirectional prepregs. This technique consists in laying and welding tapes of prepregs, building a laminate having the geometry and the desired mechanical characteristics.

This process has been widely studied [1, 2, 3, 4] with numerical models becoming more accurate, robust and complex. New techniques of computation [5, 6], based on variables' separation, have enabled models' enrichment by adding parameters as extra-coordinates [7]. Thus, the *Proper Generalised Decomposition* (PGD) method leads to multi-parametric virtual charts providing a whole set of solutions for each combination of the considered parameters [8, 9, 10, 11, 12].

Then, the computational vademecum can be exploited *on-line* for process control or process optimisation purposes. Indeed, within the AFP we want to efficiently control the heating power: tapes have to be heated enough to ensure the melting of the matrix coating the fibres and the cohesion with the previously laid tapes, while not exceeding a threshold from which material burns.

Therefore, we took advantage of the PGD to build *off-line* virtual charts in order to determine the best power associated to a draping velocity profile [13].

However, in those simulations, solutions were computed from equations provided by underlying physics of the studied phenomena, the optimisation being carried out as post-process.

We propose here to compute directly the solution of an optimisation problem in order to get the separated representations of both the field and the control to obtain it. That is to say the optimisation is made directly *off-line*,

<footnote><sup>a)</sup>n.bur@estia.fr</footnote>

<footnote><sup>b)</sup>p.joyot@estia.fr</footnote>
reducing the cost of the post-process and improving the real-time control of the AFP.

Within the next section, we present the equations governing the phenomenon under consideration. In order to have a reference solution the system is solved by standard finite element method (FEM). Thereafter, section b) focuses on the writing and solving of the optimal system. We improve the obtained results by applying Uzawa's method in section b). Lastly section b) addresses few conclusions and perspectives.

PROCESS MODELLING

AFP process can be modelled with a heat equation associated with the next boundary conditions. Heat source applies on boundary $\Gamma_P$ (see Figure 1); considering a wide enough domain $\Omega$, we set a Dirichlet's condition taking into account a fixed temperature on $\Gamma_R$; then, for the sake of simplicity, an homogeneous Neumann's condition is applied on others boundaries.

<figure><img src="image_3.png" /><figcaption>FIGURE 1. Domain of study.</figcaption></figure>

This leads to the following advection-diffusion equation 1

$$ 
\left\{
\begin{array}{ll}
\text{div}(K\nabla u) + \rho C_p \nabla u = 0 & \text{in } \Omega; \\
u = 0 & \text{on } \Gamma_R; \\
-K\partial_n u = 0 & \text{on } \Gamma_B \cup \Gamma_L \cup \Gamma_T; \\
-K\partial_n u = -\Phi & \text{on } \Gamma_P;
\end{array}
\right.
\qquad (1) $$

$$ \text{with } K = \begin{pmatrix} k_{\parallel} & 0 \\ 0 & k_{\perp} \end{pmatrix} \text{ et } V = \begin{pmatrix} v & 0 \end{pmatrix}^T. $$

This system can obviously be solved in a classical way with standard FEM. However, since we want a reference solution to latter compare with results from other methods, we have to implement comparable systems. Consequently, in order to compare solution in separated representation, we take advantage of the tensor product method introduced by R. E. Lynche, J. R. Rice and D. H. Thomas [14] [15]. The writing of discrete form of Equation (1) separating the two space coordinates and using tensor forms of operators allows the handling of shape functions in the separated representation, the system remaining solved with a global, non separated FEM.

We can thus obtain the temperature field from the input $\Phi$ modelling the laser heat flux. For some couples power/speed we compute the corresponding reference solutions.

From these results we extract the temperature profiles on boundary $\Gamma_P$ where the source is applied. Indeed, we want to control the temperature on this same boundary during the process, to reach an optimum, providing the ideal heat flux. We compute also $P_{S_{ray}}$, the power seen by $\Gamma_P$, as the integrand of the flux on this boundary multiplied by the width of the shone surface $S_{ray}$.
We gather in Table 1 some key values for these four reference solutions, for the purpose of latter comparison. Due to velocity values, we increase the number of nodes to contain the Péclet number, avoiding the requirement of a stabilisation technique

<table><caption>TABLE 1. Key values for reference solutions</caption><thead><tr><th></th><th>v = 0.001 m · s<sup>-1</sup><br>Pw = 600 W</th><th>v = 0.01 m · s<sup>-1</sup><br>Pw = 1897 W</th><th>v = 0.1 m · s<sup>-1</sup><br>Pw = 6000 W</th><th>v = 1 m · s<sup>-1</sup><br>Pw = 18 970 W</th></tr></thead><tbody><tr><th>Reference</th><td></td><td></td><td></td><td></td></tr><tr><th>U<sub>max</sub></th><td>386.076 K</td><td>371.744 K</td><td>372.283 K</td><td>372.197 K</td></tr><tr><th>P<sub>Sray</sub></th><td>320 W</td><td>1011.73 W</td><td>3200 W</td><td>10 117.3 W</td></tr></tbody></table>

## SETTING UP THE OPTIMAL SYSTEM

As announced in the Introduction, The PGD method provides multi-parametric virtual charts that can be used to control the AFP process [13].

Another way to go on is to use the optimal control theory. Within the AFP process, a heat flux provided by a laser melt the thermoplastic. The difficulty consists in determining the best power of the laser to reach an optimal temperature to melt the thermoplastic enough but not too much.

Thus we consider the following cost-function to be minimised, since the flux is applied only on $\Gamma_p$, part of the boundary

$$J(u, \Phi) = \frac{1}{2} \int_{\Gamma_p} (u - u_d)^2 + \alpha \Phi^2 d\gamma, \qquad (2)$$

subject to the state advection-diffusion equation (1), with $\alpha$ cost parameter of the command.

That way $\Phi$ is used as control and we want to reach $u_d$ on boundary $\Gamma_p$. The domain of study is depicted on Figure 1.

The corresponding Lagrangian writes

$$\mathcal{L}(u, \Phi, p) = J(u, \Phi) + \int_{\Omega} p(-\mathrm{div}(K\nabla u) + \rho C_p \nabla u) \, dx. \qquad (3)$$

To find a stationary point of the Lagrangian we set

$$\left\{ \begin{array}{l} \nabla_u \mathcal{L}(u, \Phi, p) = 0; \\ \nabla_{\Phi} \mathcal{L}(u, \Phi, p) = 0; \\ \nabla_p \mathcal{L}(u, \Phi, p) = 0. \end{array} \right. \qquad (4)$$

Expanding these equations leads to the non-linear optimality system (see [16]), whom weak form writes, with test functions $u^*$ and $p^*$,

$$\left\{ \begin{array}{l} \int_{\Omega} K \nabla u^* \cdot \nabla u + \int_{\Omega} \rho C_p \nabla u u^* - \int_{\Gamma_p} \frac{1}{a} p u^* d\gamma = 0 \quad \forall u^*; \\ \int_{\Omega} K \nabla p^* \cdot \nabla p + \int_{\Omega} \rho C_p \nabla p p^* + \int_{\Gamma_p} u p^* d\gamma = \int_{\Gamma_p} u_d p^* d\gamma \quad \forall p^*. \end{array} \right. \qquad (5)$$

In Equation (5), fields $u$ and $p$ are coupled. The computation of such a problem can be achieved using mixed formulation, within a standard FEM as well as in the PGD framework as described thereafter.

The discrete form of the state variable $u(x, z)$ and the adjoint parameter $p(x, z)$ are expressed in tensor product form

$$U = \sum_{i}^{\infty} U_{x}^{i} \otimes U_{z}^{i} \quad \text{and} \quad P = \sum_{i}^{\infty} P_{x}^{i} \otimes P_{z}^{i}. \qquad (6)$$

The vector $\Psi$ which brings together nodal values of $u$ and $p$ takes the form

$$\Psi = \begin{bmatrix} U \\ P \end{bmatrix} = \sum_{i}^{\infty} \begin{bmatrix} U_{x}^{i} \otimes U_{z}^{i} \\ P_{x}^{i} \otimes P_{z}^{i} \end{bmatrix}. \qquad (7)$$
The discretised weak form of Equation (5) is written $\Psi^* \mathcal{A} \Psi = \Psi^* \mathcal{B}$ where $\mathcal{A} = \sum_{i=1}^{8} A_x^i \otimes A_z^i$ and $\mathcal{B} = B_x \otimes B_z$ will be expressed in a tensor form, and with the test function defined by

$$ \Psi^* = \begin{bmatrix} U_x^* \otimes U_z + U_x \otimes U_z^* \\ P_x^* \otimes P_z + P_x \otimes P_z^* \end{bmatrix}. \qquad (8) $$

This discretised problem can then be solved with the PGD method, for different values of the velocity, the desired state $u_d$ coming from the corresponding reference solution.

Since the goal of optimal control is to minimise the distance between the unknown $u$ and the desired state $u_d$ on the boundary $\Gamma_p$, we retrieve the wanted temperature profile on this boundary. However, to reach it, the flux, computed by PGD, to be applied take substantially away compared to the FEM input, as speed increases (see Figure 2).

<figure><img src="image_4.png" /></figure>

$v = 0.001 \text{ ms}^{-1}$

<figure><img src="image_5.png" /></figure>

$v = 0.01 \text{ ms}^{-1}$

<figure><img src="image_6.png" /></figure>

$v = 0.1 \text{ ms}^{-1}$

<figure><img src="image_7.png" /></figure>

$v = 1 \text{ ms}^{-1}$

<caption>FIGURE 2. Flux on $\Gamma_p$ for different speeds.</caption>

## UN-MIXING THE OPTIMALITY SYSTEM

Applying the PGD framework on the optimality system (5) can produce, in our case, weird results. Instead of searching a workaround within numerical stabilisation techniques, we propose to solve this system without mixed elements. For this purpose, we use Uzawa technique [17] [18].

Given an initial guess $p^{(0)}$ for $p$, Uzawa's method consists in our case of the following coupled iteration:

$$ \begin{cases} \mathcal{A}^u u^{(k+1)} = \mathcal{A}^{pu} p^{(k)} \\ p^{(k+1)} = p^{(k)} + \omega (\mathcal{A}^{pu} (u^{(k+1)} - u_d) + \mathcal{A}^p p^{(k)}) \end{cases} \qquad (9) $$
where $\omega > 0$ is a relaxation parameter. Tensors are defined by

$$A^u = k_{\parallel} K_x^u \otimes M_z^u + M_x^u \otimes k_{\perp} K_z^u + \rho C_p \nu H_x^u \otimes M_z^u$$

$$A^p = k_{\parallel} K_x^p \otimes M_z^p + M_x^p \otimes k_{\perp} K_z^p - \rho C_p \nu H_x^p \otimes M_z^p$$

$$A^{pu} = -\frac{1}{\alpha} M_x^{pu} \otimes M_z^{pu}$$

$$A^{up} = M_x^{up} \otimes M_z^{up}$$

This coupled iteration is computed within a fixed-point loop, each field u and p being solved by PGD. As previously, Figure 3 shows the computed flux to reach the desired temperature $u_d$ on the boundary $\Gamma_p$.

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 3. Flux on $\Gamma_p$ with Uzawa's method</figcaption></figure>

Table 2 summarises the previous simulations, collecting some key-values. Thus, solving the optimality system with the PGD algorithm, but without mixed elements using Uzawa method leads to consistent results, since both temperature and heat flux remain similar to those expected (using standard FEM).

## CONCLUSION

In order to control the AFP process, we proposed a method based on the PGD technique allowing the inclusion of parameters as extra-coordinates. Due to numerical instabilities, we transformed a mixed formulation to a coupled problem, improving significantly the results.

Next step consists in considering the desired state $u_d$ as a new coordinate. Thus we will be able to build directly the process control as a multi-parametric virtual chart.

## REFERENCES

[1] F. O. Sonmez, H. T. Hahn, and M. Akbulut, Journal of Thermoplastic Composite Materials 15, 525–544 (2002).
<table><caption>TABLE 2. Summary of key values</caption><thead><tr><th></th><th>ν = 0.001 m · s<sup>-1</sup><br>Pw = 600 W</th><th>ν = 0.01 m · s<sup>-1</sup><br>Pw = 1897 W</th><th>ν = 0.1 m · s<sup>-1</sup><br>Pw = 6000 W</th><th>ν = 1 m · s<sup>-1</sup><br>Pw = 18 970 W</th></tr></thead><tbody><tr><td colspan="5"><strong>Reference - kFEM</strong></td></tr><tr><td>U<sub>max</sub></td><td>386.076 K</td><td>371.744 K</td><td>372.283 K</td><td>372.197 K</td></tr><tr><td>P<sub>S,ray</sub></td><td>320 W</td><td>1011.73 W</td><td>3200 W</td><td>10 117.3 W</td></tr><tr><td>PGD</td><td>16 modes</td><td>183 modes</td><td>200 modes</td><td>200 modes</td></tr><tr><td>U<sub>max</sub></td><td>386.077 K</td><td>384.143 K</td><td>364.924 K</td><td>367.761 K</td></tr><tr><td colspan="5"><strong>Ctrl - kFEM</strong></td></tr><tr><td>U<sub>max</sub></td><td>385.542 K</td><td>374.386 K</td><td>372.249 K</td><td>372.164 K</td></tr><tr><td>P<sub>S,ray</sub></td><td>319.919 W</td><td>1011.69 W</td><td>3199.99 W</td><td>10 117.3 W</td></tr><tr><td>Ctrl - PGD</td><td>26 modes</td><td>156 modes</td><td>200 modes</td><td>200 modes</td></tr><tr><td>U<sub>max</sub></td><td>387.108 K</td><td>371.275 K</td><td>369.883 K</td><td>394.982 K</td></tr><tr><td>P<sub>S,ray</sub></td><td>337.083 W</td><td>-4177.73 W</td><td>119 749 W</td><td>1.650 41 × 10<sup>6</sup> W</td></tr><tr><td colspan="5"><strong>Ctrl - Uzawa PGD</strong></td></tr><tr><td>U</td><td>17 modes</td><td>80 modes</td><td>80 modes</td><td>80 modes</td></tr><tr><td>U<sub>max</sub></td><td>383.359 K</td><td>381.267 K</td><td>400.239 K</td><td>360.067 K</td></tr><tr><td>P</td><td>6 modes</td><td>5 modes</td><td>3 modes</td><td>3 modes</td></tr><tr><td>P<sub>S,ray</sub></td><td>318.585 W</td><td>1015.02 W</td><td>3238.99 W</td><td>10 087.2 W</td></tr></tbody></table>

[2] R. Pitchumani, S. Ranganathan, R. C. Don, J. W. Gillespie, and M. A. Lamontia, International Journal of Heat and Mass Transfer **39**, 1883–1897 (1996).

[3] R. Schledjewski and M. Latrille, Composites Science and Technology **63**, 2111–2118 (nov 2003).

[4] M. A. Lamontia, M. B. Gruber, J. Tierney, J. W. Gillespie, B. J. Jensen, and R. J. Cano, Jec-Sampe 1–8 (2009).

[5] A. Ammar, B. Mokdad, F. Chinesta, and R. Keunings, Journal of Non-Newtonian Fluid Mechanics **139**, 153–176 (2006).

[6] A. Ammar, B. Mokdad, F. Chinesta, and R. Keunings, Journal of Non-Newtonian Fluid Mechanics **144**, 98–121 (2007).

[7] F. Chinesta, A. Leygue, B. Bognet, C. Ghnatios, F. Poulhaon, F. Bordeu, A. Barasinski, A. Poitou, S. Chatel, and S. Maison-Le-Poec, International Journal of Material Forming **7**, 81–92 (sep 2014).

[8] É. Prulière, F. Chinesta, and A. Ammar, Mathematics and Computers in Simulation **81**, 791–810 (2010).

[9] B. Bognet, F. Bordeu, F. Chinesta, A. Leygue, and A. Poitou, Computer Methods in Applied Mechanics and Engineering **201-204**, 1–12 (2012).

[10] C. Ghnatios, F. Chinesta, E. Cueto, A. Leygue, A. Poitou, P. Breitkopf, and P. Villon, Composites Part A: Applied Science and Manufacturing **42**, 1169–1178 (may 2011).

[11] C. Ghnatios, F. Masson, A. Huerta, A. Leygue, E. Cueto, and F. Chinesta, Computer Methods in Applied Mechanics and Engineering **213-216**, 29–41 (mar 2012).

[12] F. Chinesta, A. Leygue, F. Bordeu, J. V. Aguado, E. Cueto, D. González, I. Alfaro, A. Ammar, and A. Huerta, Archives of Computational Methods in Engineering **20**, 31–59 (jan 2013).

[13] N. Bur, P. Joyot, C. Ghnatios, P. Villon, E. Cueto, and F. Chinesta, Advanced Modeling and Simulation in Engineering Sciences **3**, p. 4 (dec 2016).

[14] R. E. Lynch, J. R. Rice, and D. H. Thomas, Bulletin of the American Mathematical Society **70**, 378–384 (1964).

[15] R. E. Lynch, J. R. Rice, and D. H. Thomas, Numerische Mathematik **6**, 185–199 (dec 1964).

[16] J.-L. Lions, *Optimal Control of Systems Governed by Partial Differential Equations* (Springer-Verlag, 1971) p. 396.

[17] H. Uzawa, in *Studies in Linear and Nonlinear Programming*, edited by K. J. Arrow, L. Hurwicz, and H. Uzawa (Stanford University Press, Stanford, California, 1958), pp. 154–165.

[18] R. Simon, “Multigrid Solvers for Saddle Point Problems in PDE-Constrained Optimization Dissertation,” Ph.D. thesis, Johannes Kepler Universität Linz 2008.