Michel Lemaire
email: michel.lemaire1@uqtr.ca

Daniel Massicotte
email: daniel.massicote@uqtr.ca

Jeremy Poupart
email: jeremy.poupart@uqtr.ca

Patrice Quinton
email: patrice.quinton@ens-rennes.fr

Sanjay Rajopadhye
email: sanjay.rajopadhye@colostate.edu

Polyhedra at Work: Automatic Generation of VHDL Code for the Sherman-Morrison Formula

Keywords: Simulation of Electrical Circuits, Polyhedral Model, Sherman-Morrison Formula, vhdl, fpga

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The simulation of electrical circuits in real-time is used in many applications [START_REF] Belanger | The what, where and when of´real-time simulations[END_REF][START_REF] Farid | Circuit Simulation[END_REF]. Real-time simulations are often needed, for example, when some parts of the system under design are replaced by a simulator.

In most methods, electrical simulation amounts to solving a linear system 𝐴𝑥 = 𝑏 where 𝐴 is the admittance matrix of the circuit, 𝑏 is the values of currents and voltages of the circuit at a given instant of time, and 𝑥 is (the vector of) the current and voltages after a small time step. The values of the 𝐴 matrix represent the admittance (reciprocal of the resistances, basically) of the edges of the circuit, represented as a graph. One step of simulation is to compute 𝑥 = 𝐴 -1 𝑏, which requires to compute the inverse of 𝐴.

However, electrical circuits contain switches, the admittance of which may be modified from one step to the next one, if the status (open vs. closed) of the switch changes. Thus, matrix 𝐴 -1 may change from one step to the next.

IMPACT 2024, January 17, 2024, Munich, Germany 2024.

For real-time simulations, current electrical simulators are implemented on fpga platforms which provide the low latency needed as well as resources to execute in parallel some calculations. In practice, when the size of the circuit is small enough, inverse matrices 𝐴 -1 are pre-computed and stored on the memory of the fpga, but this is of course limited by the size of the circuit, and by the number of switches that it contains.

An alternative is to update the inverse of the 𝐴 matrix when the state of some switches changes. But, without optimization, the 𝑂 (𝑁 3 ) complexity of matrix inversion prevents such a method to be done in real-time, on the current fpga platforms.

Implementing such an algorithm on a fpga is not an easy task, even with modern high-level synthesis tools that are more an more often used. This task requires the exploration of several solutions, to make sure that under the limited resources of the fpga platform, the performance of the design meets the requirements of the application. Besides, the numerical stability of the chosen algorithm has to be ensured, since simulations may last several hours or even days.

In this paper, we explore the automatic generation of synthesizable vhdl code for the 𝑂 (𝑁 2 ) Sherman-Morrison algorithm [START_REF] Sherman | Adjustment of an Inverse Matrix Corresponding to a Change in One Element of a Given Matrix[END_REF] for incremental matrix inversion, in the context of electrical circuit simulation, using an approach based on the Polyhedral Equational Model. In this model, calculations are expressed as recurrence equations that can be seen as single-assignment statements describing the problem to solve. We use tools based on the alpha language [START_REF] Verge | The ALPHA Language and its Use for the Design of Systolic Arrays[END_REF][START_REF] Rajopadhye | Alphabets: An Extended Polyhedral Equational Language[END_REF], one of the long-standing vehicles for this approach. Our contribution is to present a fully automatic design flow that starts from an alpha description and produces in a few seconds, a vhdl program which can be synthesized for a fpga.

This paper is organized as follows. In Section 2, we present the context of this research, the simulation of electrical circuits. Section 3 describes the Sherman-Morrison formula. The Polyhedral Equational Model is presented in Section 4, together with the alpha language and associated tools. Section 5 describes a first implementation of the Sherman-Morri-son algorithm using alpha. The translation to vhdl is explained in Section 6. In Section 7, we present an optimized version of the Sherman-Morrison formula, in the context of the electrical simulation. Section 8 is devoted to results. In Section 10, we discuss some aspects of this work, and Section 11 concludes and presents some future research directions.

Simulation of Electrical Circuits

Fig. 1 shows the schematics of the electrical circuit of a power converter.

Simulating the behaviour of such a circuit amounts to solve the linear system 𝐴𝑥 = 𝑏 to obtain, from the current voltage and current vector 𝑏, the new values 𝑥 of these quantities.

When the electrical circuit contains switches-as is the case in this power-converter,-the 𝐴 matrix may change between two simulation steps, since switches are modeled by a resistance, very high when the circuit is open, and very low otherwise. The matrix 𝐴, being dependant on the switches states, requires an inversion at every switch state change.

The real-time simulation of such a circuit can be done on a fast accelerator, for example, a fpga platform, but to reach the low latency required by real-time, one has either to store all possible values of the 𝐴 -1 , which can be done if the circuit is small enough, or to re-compute 𝐴 -1 as needed.

In general, the complexity of solving a linear system of size 𝑁 is 𝑂 (𝑁 3 ). However, in the particular case of a circuit simulation, only a few elements of the 𝐴 matrix change.

In the example shown by Fig. 2, the 𝑔 1 value corresponds to the admittance of switch 𝑆 1 of Fig. 1, and this value only appears in elements 𝑎 11 , 𝑎 31 , 𝑎 13 and 𝑎 33 of matrix 𝐴.

Therefore, it is interesting to consider methods to update the 𝐴 matrix with a lower complexity, among which, the Sherman-Morrison approach.

The Sherman-Morrison Formula

Given a square, invertible matrix 𝐴, the Sherman-Morrison algorithm allows one to efficiently compute the inverse of an order-one perturbation of 𝐴, according to the formula:

(𝐴 + 𝑢𝑣 𝑇 ) -1 = 𝐴 -1 -𝜎𝐴 -1 𝑢𝑣 𝑇 𝐴 -1 , ( 1
)
where

𝜎 = 1 1 + 𝑣 𝑇 𝐴 -1 𝑢 . ( 2
)
In this formula, the outer-product 𝑢𝑣 𝑇 represents an orderone perturbation of the matrix 𝐴.

Computing this formula requires the calculation of two vector-products-𝐴 -1 𝑢 and 𝑣 𝑇 𝐴 -1 -that require each 𝑂 (𝑁 2 ) operations, an outer product, and a subtraction of matriceseach one also 𝑂 (𝑁 2 ) operations-plus a dot product and a division. This is better than a full-matrix inversion which needs 𝑂 (𝑁 3 ) operations.

Note that higher-order perturbations of 𝐴 can be dealt with by successive applications of the Sherman-Morrison method.

Several problems have to be faced, when implementing this algorithm in the context of electrical simulation, in practice:

• The latency of the algorithm, i.e., the number of steps between the moment that the changes of 𝐴 are known, and the moment when the update of 𝐴 -1 is computed, must be very small, which requires a highly parallel implementation and the use of a special-purpose platform such as a FPGA. • The accuracy of the implementation must be such that the simulation results remain correct, even after long simulations.

• Many optimisations might have to be tried, necessitating rapid design-space exploration.

The Polyhedral Equational Model

The Polyhedral Equational Model is a formalism based on recurrence equations which has been developed over years to represent regular calculations and transform these representations in order to derive executable code for parallel architectures, or even, to produce hardware descriptions [START_REF] Quinton | Automatic Synthesis of Systolic Arrays from Uniform Recurrent Equations[END_REF][START_REF] Sanjay | Synthesizing systolic arrays from recurrence equations[END_REF].

The alpha language [START_REF] Verge | The ALPHA Language and its Use for the Design of Systolic Arrays[END_REF][START_REF] Mauras | Alpha : un langage équationnel pour la conception et la programmation d'architectures parallèles synchrones[END_REF][START_REF] Rajopadhye | Alphabets: An Extended Polyhedral Equational Language[END_REF] that we shall consider in this paper is a functional language invented to express polyhedral equations. In this section, we first present briefly the alpha language (4.1), then its transformations (4.2), and the tools which are available (4.3).

The alpha language

Fig. 3 shows the alpha description of a matrix-vector algorithm that will be used in the following. After the key-word system, the name of the program is given, followed by the definition of a size parameter 𝑁 with its domain.

The inputs of the program (lines 2 and 3) are a matrix 𝑎 and a vector 𝑣 (for the sake of simplicity, of integers). The output (line 4) is a vector 𝑐. The domains of these variables are given as polyhedra, in the obvious way.

The equations that define the calculations make use of an intermediate-local-variable 𝑋 that is defined using a simple induction (line 8 to 12), and the definition of the output 𝑐 (line 13) as the last element of this recurrence.

This program is shown in its so-called array form, that allows classical numerical calculations to be expressed in a conventional way. However, the definition of alpha is functional, and the basic constructs of the language allow formal manipulations of code suitable for parallel expression to be performed [START_REF] Mauras | Alpha : un langage équationnel pour la conception et la programmation d'architectures parallèles synchrones[END_REF].

The alpha language contains the notion of subsystems that allows a system to be reused in a hierarchy of definitions, while maintaining the properties of the core language.

For example, the matrix-vector alpha model shown in Fig. 3 can be included inside another system by:

use matVect[N] (A, x) returns (Y)
provided that the inputs 𝐴, 𝑥 and the output 𝑌 are declared with the same domains and scalar types as the formal parameters 𝐴, 𝑣 and 𝑐 in the matVect model. The matvect[N] part of the statement indicates that the model is used with the same parameter 𝑁 . In general [START_REF] Florent De Dinechin | Structuration of the ALPHA language[END_REF], the form of such a statement includes a mapping domain, where the model can be instanciated, and the parameter of the instanciated system can be any affine expression of the size parameters of the calling system.

Main Polyhedral Transformations

The following semantics-preserving transformations of code will be extensively used during the synthesis:

• Substitution aka., inlining: a variable can be replaced by its definition, anywhere, by simple cut and paste. Conversely, by anti-substitution, one can replace any expression by a new variable, whose definition domain can be computed easily from the constituents of the expression. • Normalization: any expression can be rewritten in the simple case-restriction-dependence form, much as shown in the definition of 𝑋 in the program of Fig. 3.

• Change-of-basis: the domain of definition of a local

variable can be replaced by its image under a linear, unimodular mapping, by a simple, syntactic rewriting. This allows space-time mapping to be done: provided the mapping (𝑧 → 𝑡, 𝑝) that transforms the indexes 𝑧 of some variable 𝑉 into a pair of time and processor coordinates, the new program becomes executable, by evaluating a calculation 𝑉 (𝑧) on processing unit 𝑝 at time 𝑡. • Pipelining: by identification of values that are common to several expressions, reuse of these common values can be performed through the rephrasing of some equations using linear inductions. • Simplification of reductions: the alpha language allows the direct expression of reductions such as Σ, Π, etc. For example, we could have written the definition of 𝑐 in the matrix-vector program as

c[i] = reduce( +, j, a[i,j] * v[k]) that reads obviously 𝑐 [𝑖] = 𝑗 𝑎[𝑖, 𝑗] × 𝑣 [ 𝑗].
• Inlining and exlining: the use of a model can be replaced by its core definition, provided some simple composition rules are applied to the equations. This inlining process is most often used prior to the transformation of a design. However, a fully structured [START_REF] Quinton | Structured Scheduling of Recurrence Equations: Theory and Practice[END_REF] design path allows a design to be done in a hierarchical way, but this will not be used in this paper examples. Conversely, exlining allows some calculations to be separated from a given system, and this operation is actually used in the mmalpha design trajectory to separate the control and datapath parts of the program, since they are translated in a different way in vhdl.

Equational Polyhedral Tools

There exists currently two sets of tools allowing alpha programs to be transformed. Both are based on the same abstract syntax, although, for historical reasons, their concrete syntax slightly differs. mmalpha [START_REF] Charot | Modeling and Scheduling Parallel Data Flow Systems using Structured Systems of Recurrence Equations[END_REF] was developed at Irisa over years, and it targets, currently, the generation of fine-grain architectures. It includes the following transformations:

1. Inlining of sub-systems. 2. Scheduling (using affine by variable schedules). 3. Apply a schedule to operate a time-space mapping of the program. 4. Find pipelines. 5. Rewrite a scheduled and allocated program into a RTLlike code suitable for hardware generation. 6. Generation of naive, single-assignment C code, for simulation and verification purpose. 7. Generation of synthesizable vhdl code. 8. Generation of stimuli files for the hardware description and simulation.

The alphaz tool [START_REF] Yuki | AlphaZ: A System for Design Space Exploration in the Polyhedral Model[END_REF] is merely a transformation exploration system, where the user can provide a script that allows code to be generated for a larger class of parallel architectures. An important difference with mmalpha is that alphaz allows:

1. Tiling, i.e., organizing the calculations in blocks of the same size called tiles, where the executions follows a sequential order. This is essential to obtain good tradeoffs between time and space, in order to optimize the resource utilization.

2. Optimized code generation, including memory mapping of variables [START_REF] Bastoul | Code Generation in the Polyhedral Model Is Easier Than You Think[END_REF][START_REF] Quilleré | Optimizing Memory Usage in the Polyhedral Model[END_REF]. 3. Simplification of the reductions. In high dimensional designs, it may happen that sharing the calculation of sub-expressions allows the complexity of the whole program to be reduced by several orders of magnitude [START_REF] Gautam | Simplifying Reductions[END_REF].

It is possible to combine both mmalpha and alphaz in order to organize the design path of an application, since at any time, the design can be expressed as an abstract alpha program combined with some annotations expressing the current state of transformations.

Expression and Synthesis of the Simulation Algorithm

In this section, we present the alpha code for an implementation of the Sherman-Morrison formula (5.1), then we describe the synthesis flow of this program (5.2).

alpha Code

A naive description of the Sherman-Morrison algorithm using the alpha language is shown in Fig. 4. This encoding maps almost directly the equations of the definition into the language. Lines 1 to 4 are include statements that allow some other, basic functions to be used. The outProd function computes the outer products of two vectors 𝑥 and 𝑦 of size 𝑁 into the square matrix (𝑥 𝑖 𝑦 𝑗 ) 1≤𝑖,𝑗 ≤𝑁 .

The inputs are the matrix 𝐵 that contains the current version of 𝐴 -1 , and the vectors 𝑢 and 𝑣. The output newB is the updated matrix.

Remark:

In this version, as well as throughout this paper, the 𝜎 coefficient is set to 1, in order to avoid a division. Indeed, the use of a division can be avoided by a precalculation of all its possible values, and then, memorized in a look-up table. This part of the design is not considered here, in order to simplify the presentation.

The Synthesis Process

In this section, we illustrate the use of mmalpha to generate vhdl code, for various situations resulting from several optimization choices.

The goal of the synthesis is to derive a vhdl block that can be eventually included in a larger schematics, targeted to the Mathworks System Generator tool [START_REF]Xilinx 2022. System Generator for DSP Overview (Online)[END_REF].

Principles of the Synthesis.

The synthesis process is fully automatic. It consists of a series of well-identified steps which require no a priori information others than the alpha program.

Until the last, vhdl generation step, the synthesis results in a transformation of the initial alpha code that is semantically equivalent to the initial description. 1. The initial program is parsed, and calls to subsystems are inlined. Another possibility is to do the synthesis in a structured way, without inlining subsystems, but this facility is not used in this example. 2. The parsing includes a thorough type-checking of the alpha code, in particular, the verification that domains of expressions are consistent, and that calls to subsystems are coherent. 3. The alpha program is then scheduled, using the socalled vertex method [START_REF] Mauras | Scheduling Affine Parameterized Recurrences by Means of Variable Dependent Timing Functions[END_REF][START_REF] Sanjay | Synthesizing systolic arrays from recurrence equations[END_REF][START_REF] Rajopadhye | On Synthesizing Systolic Arrays from Recurrence Equations with Linear Dependencies[END_REF] In addition, the program is separated into three parts: a wrapper, having the same signature that the initial alpha program, a hardware module comprising the data-path equations, and a controller, containing the equation that define boolean control signals. This representation is again, semantically equivalent to the initial alpha program. 6. All transformations are done independently on the size parameters. Before the generation of the vhdl code, the parameters are set to the value required for the design. 7. Finally, the vhdl code is generated. We detail this last step in section 6.

Simulations.

The synthesis process contains several ways to simulate a design under synthesis, at various levels. These simulations are intended to help verifying the initial alpha program, to generate stimuli for the vhdl code, and also, to test the mmalpha tool. A first simulation is based on a strict application of the denotational semantics of alpha, as defined by Mauras [START_REF] Mauras | Alpha : un langage équationnel pour la conception et la programmation d'architectures parallèles synchrones[END_REF] and De Dinechin [START_REF] Florent De Dinechin | Structuration of the ALPHA language[END_REF]. It allows an early verification of the design, as well as generation of random inputs for creating automatically stimuli files.

Another simulation is done by a translation of the alpha program into a C program using Quillere's translator [START_REF] Quilleré | Optimizing Memory Usage in the Polyhedral Model[END_REF].

Finally, another simulation is based on repeated substitutions and normalization of equations. It is very slow, but can be applied to the final output of the synthesis, and is very helpful to verify the mmalpha tool.

Translation to vhdl

The translation to vhdl targets a simple synthesizable form. At the end of the synthesis, the alpha program is separated into a control subsystem, a data-path subsystem, and a wrapper subsystem. The control subsystem is translated into a vhdl finite state machine. The wrapper subsystem is used to describe the detailed input/output of the system, and to generate a simulation script. In this section, we concentrate on the generation of the vhdl code for the data-path part. In 6.1, we explain how alpha variables are translated into vhdl arrays, and in 6.2, we describe the translation of the data-path equations. It should be noticed that the translation process that we describe here does not consider the input/output and the memorization of data. We report the discussion of this important limitation to Section 11.

Array Types of Variables

Variables are mapped to vhdl arrays. Any variable of the hardware module has a domain that is time-space mapped, of the form {𝑡, 𝑝 | 𝑐} where 𝑐 is a set of linear constraints. The bounding box of such a domain can be easily found, and it is used to define the vhdl array type of this variable.

This process is sufficient for most of the variables, provided their life-time is bounded.

The notion of life-time was introduced in Feautrier and Lefebvre [START_REF] Lefebvre | Automatic storage management for parallel programs[END_REF], and Quilleré and Rajopahdye [START_REF] Quilleré | Optimizing Memory Usage in the Polyhedral Model[END_REF], to cope with the memory-mapping of variables before the generation of parallel code.

The life-time of a variable instance, say 𝑋 [𝑡, 𝑝] is the longest amount of time that this instance is required to be memorized to be consumed by the calculation of any other instance 𝑌 [𝑡 ′ , 𝑝 ′ ].

Finding out the life-time of variables amounts to compute the maximum value of 𝑡 ′ -𝑡 over the set of all these points, which can be done thanks to parameter integer programming, for example, using the Pip software [START_REF] Feautrier | Parametric integer programming[END_REF]. This maximum, in general, depends on 𝑝, and on the size parameters (here, 𝑁 ). However, as we do the generation of vhdl code after setting the value of the parameters, and since 𝑝 is bounded, the life-time is a constant 1 .

Several situations may occur.

• If the life-time is 0, the instance 𝑋 [𝑡, 𝑝] is always used at the same instant of time when it is produced, and no memorization is required. • If the life-time is 1, we implement this as a simple register. • If the life-time is strictly higher than 1, we implement the memorisation of 𝑋 by adding an additional variable, implementing a registered version of 𝑋 , with an extra dimension to the declaration array of 𝑋 , the size of which is defined by the life-time value.

Translation of the equations

The equation defining a variable is unique, and has the form X[t,p] = ..., where 𝑝 is a set of indexes representing the space dimensions of 𝑋 . This equation is translated into a nest of vhdl For generate statements, whose bounds are computed from the domain of 𝑝 thanks to a domain to loop computation algorithm [START_REF] Ancourt | Scanning Polyhedra with DO Loops[END_REF][START_REF] Quilleré | Optimizing Memory Usage in the Polyhedral Model[END_REF]. The body of this nest depends on the type of equation. For a combinatorial equation, the translation is almost direct, by mapping the alpha operators to those of vhdl.

For a multiplexer, the translation is also very simple, adding a When condition to represent the condition.

The translation of a simple connection is also very simple. More difficult is the translation of so-called register equations, which raise actually both memorization and communication issues.

Consider such an equation

X[t,p] = Y[f(t,p),g(t,p)]
where 𝑓 (𝑡, 𝑝) is the expression of the time in 𝑌 and 𝑔(𝑡, 𝑝) is the expression of the space. Here, 𝑡 and 𝑝 are vectors of indexes, and 𝑓 and 𝑔 are multi-dimensional affine functions.

When 𝑓 is the identity, such an equation represents a communication between the processors of 𝑋 and those of 𝑌 . In the current implementation, we translate directly the 𝑔 function in vhdl code, letting the vhdl hardware synthesis tools cope with the communication implementation. In some extreme situations, this may result in bad hardware implementations, and further optimizations based on the analysis of polyhedral equations is required (see for example Van Dongen and Quinton [START_REF] Quinton | The Mapping of Linear Recurrence Equations on Regular Arrays[END_REF]).

When 𝑓 is not the identity, we also separate the situations where 𝑓 (𝑡, 𝑝) depends or not on 𝑝. In the most general case, we compute the maximum value of the difference 𝑡 -𝑓 (𝑡, 𝑝), say 𝑚 and we implement the access to the proper value 𝑌 [𝑓 (𝑡, 𝑝), 𝑔(𝑡, 𝑝)] through the use of a shift-register that stores the 𝑚 most recent values of 𝑌 .

Although this implementation may not be optimal, it is guaranteed to produce correct code.

Optimized Code for the Electrical Simulation

Although the description of the Sherman-Morrison formula of Fig. 3 could be used to solve the updating of the admittance matrix, some simplifications can be made to this program by taking into account the structure of this problem.

The 𝐴 matrix elements that are changed, when the status of a switch changes, are the 4 elements 𝑎 𝑖𝑖 , 𝑎 𝑖 𝑗 , 𝑎 𝑗𝑖 and 𝑎 𝑗 𝑗 , and the corresponding modification is given by the matrix

0 0 0 0 0 0 -𝑑 𝑖 0 𝑑 𝑗 0 0 0 0 0 0 0 𝑑 𝑗 0 -𝑑 𝑖 0 0 0 0 0 0 , ( 3
)
which can be represented by the product 𝑢𝑣 𝑇 where 𝑢 = (0, 𝑑 𝑖 , 0, -𝑑 𝑗 , 0) and 𝑣 = (0, -1, 0, 1, 0). The modification of the 𝐴 -1 matrix, given by

𝐴 -1 ← 𝐴 -1 + 𝜎𝐴 -1 𝑢𝑣 𝑇 𝐴 -1 , (4)
can be done in the following way:

1. Select the columns 𝑐 𝑖 and 𝑐 𝑗 of 𝐴 -1 , corresponding to the numbers 𝑖 and 𝑗 of the non-zero elements of 𝑢 and 𝑣.

The product 𝑙

= 𝐴 -1 𝑢 is then equal to 𝑑 𝑖 × 𝑐 𝑖 -𝑑 𝑗 × 𝑐 𝑗 . 3. The product 𝑟 = 𝑣 𝑇 𝐴 -1 is equal to 𝑐 𝑖 + 𝑐 𝑗 . 4. 𝐴 -1 𝑢𝑣 𝑇 𝐴 -1 is then 𝜎 (𝐴 -1 -𝑙𝑟 𝑇 ).

Results

We have generated the vhdl code for two versions of the Sherman-Morrison formula, first, the full version of Fig. 3, then the optimized version of Section 7 2 .

The produced vhdl code was simulated and then transformed into a schematics using the Vivado System Generator software [START_REF]Xilinx 2022. System Generator for DSP Overview (Online)[END_REF]. Table 1 displays the results of these experiments.

In this table, the number of processors represents the maximum space size of the alpha equations, at the end of the design, before the vhdl translation. However, the actual vhdl design is not organized as a collection of identical processing units, as each equation produces an independent generate statement. As a consequence, the number of DSPs of the design depends on resource usage optimization done by the fpga tools.

The full Sherman-Morrison algorithm3 produces a spacelinear design, with linear latency. The complexity is dominated by the matrix-vector product, which requires 𝑁 processing units.

The optimized algorithm produces a design with 𝑁 2 processing units required to update the 𝐴 -1 matrix, which is kept in registers. The latency is constant, independent on the size of the design.

As expected, the synthesis time-i.e., the time needed to produce the vhdl code-is independent on the size parameter 𝑁 .

The vhdl code, generated directly from the final alpha equations, is very readable, as each statement is a direct translation of an alpha equation, as shown by the example of Appendix B.

Other Works

Implementations of the Sherman-Morrison formula to accelerate Electrical simulation on fpga has been considered by several authors (see [START_REF] Hadizadeh | A Matrix-Inversion Technique for FPGA-based Real-Time EMT Simulation of Power Converters[END_REF] for example). The systematic exploration and optimization of fpga architectures to implement this method is more recent, and either is based on handwritten vhdl code as in [START_REF] Hadizadeh | A Matrix-Inversion Technique for FPGA-based Real-Time EMT Simulation of Power Converters[END_REF], or on C code annotated for the use of high-level-synthesis tools, as in [START_REF] Poupart | FPGA Implementation of Sherman-Morrison Formula Using High-Level Synthesis and Graphical Blocks Programming[END_REF].

The generation of architectures from Polyhedral models has been less addressed than the generation of parallel programs. Most research on the polyhedral model start from affine loops in imperative programs, and not on equations. Most often, the target is the production of parallel programs, but the success of high-level synthesis tools have led some projects to target the generation of C code adapted to fpga compilation [START_REF] Devos | Finding and Applying Loop Transformations for Generating Optimized FPGA Implementations[END_REF][START_REF] Zhao | POLSCA: Polyhedral High-Level Synthesis with Compiler Transformations[END_REF][START_REF] Zuo | Improving polyhedral code generation for highlevel synthesis[END_REF][START_REF] Zuo | Improving High Level Synthesis Optimization Opportunity through Polyhedral Transformations[END_REF].

Recent research on parallel accelerators for machine learning consider systolic arrays as a target [START_REF] Wang | Search for Optimal Systolic Arrays: A Comprehensive Automated Exploration Framework and Lessons Learned[END_REF]), and make use of polyhedral optimization to obtain efficient designs.

Research on the polyhedral equational approach to target vhdl has been pursued in [START_REF] Charot | Modeling and Scheduling Parallel Data Flow Systems using Structured Systems of Recurrence Equations[END_REF][START_REF] Quinton | Efficient hardware implementation of data-flow parallel embedded systems[END_REF], and the research presented in this paper results of a large extension of the mmalpha software done to implement such methods.

Discussion

Several aspects of this research deserve comments.

• The choice of implementing the computation of the 𝜎 coefficient using a look-up table is motivated by the difficulty of computing the division on a fpga, and the observation that the number of values taken by this coefficient is very limited (a few hundreds of values, depending on the size of the 𝐴 matrix). However, this depends largely on the fpga platform that is used for the implementation, and such a choice might be reconsidered if needed. 1. Results for the synthesis of several versions of the Sherman-Morrison formula, showing the size of the 𝐴 matrix, the fpga resources used by the design (number of DSPs, of Look-up-tables, and of flip-flops), the latency of the design in number of cycles, the number of processors of the design, the number of lines of alpha code and the number of vhdl lines of the code generated.

• One may worry about the numerical stability of the Sherman-Morrison method used here. This issue, shared by all matrix inversion methods, can be solved by several approaches. In the context of an fpga implementation, one can try to find out, experimentally, the best number representation, either as fixed-point or floating-point. Another flexibility is given by the possibility to choose the initial 𝐴 -1 matrix among a set of carefully selected ones, in such a way that the error of the inversion is minimized. All these aspects are currently under investigation. • The data types of the alpha variables can be chosen according to the needs of the algorithm being implemented: integers, fixed point or floating point numbers of various size. In fact, the whole synthesis process is largely independent on the choice of these datatypes, which matters only when generating the vhdl code. • The synthesis method that is described here has been applied successfully to many other algorithms: matrix operations, string matching, signal processing algorithms for example. It has currently been validated by a functional simulation of the vhdl code. • Due to lack of space, we do not detail here, the time needed by each step of the synthesis. The scheduling step, which is supposed to be the critical step of the synthesis, takes only a small fraction of the total synthesis time, on the order of 10%, using the linear solver of the mathematica software. Most of the time is spent in the transformation of the time-space mapped alpha code into a multi-dimensional RTL description, which is essentially due to the use of the interpreted Wolfram language. • In its current status, the mmalpha software does not provide a facility to generate input/outputs and memorization of data. In the context of the present application, this is not a limitation, since the generated code is meant to be included in a schematic description used by the System Generator tool. The solution to overcome this limitation would be to implement transformations allowing I/O and memorization to be performed. In general, obtaining efficient implementations is a difficult problem, as it depends heavily on the hardware platform that is available.

Conclusion and Future Work

We have shown how the Sherman-Morrison formula for recomputing the inverse of a square matrix subject to a oneorder perturbation can be expressed using a Polyhedral equational description, and then, transformed into a vhdl hardware description suitable for fpga implementation.

We have presented results for the full Sherman-Morrison algorithm, and for a version optimized for the simulation of electrical circuits.

Our results show that the synthesis time is independent of the size of the problem, which is a main advantage of the polyhedral model, where calculations, be they expressed as loops or as multi-dimensional equations, are represented by operations on collections of data.

We have shown that a representation by polyhedral equations allows efficient parallel designs to be obtained automatically.

Future directions of research are to enlarge the set of available transformations, in order to target various architecture organizations, on fpga on or other types of parallel architectures. To this end, an interconnection between mmalpha and alphaz would represent an interesting experimentation platform.

Another direction would be to take advantage of the functional nature of the alpha languages to explore algebraic transformations-for example, linear algebra properties-since the rewriting of equations certainly offers large opportunities of code simplifications.

A Implementation of mmalpha mmalpha is implemented in the Wolfram language, under the mathematica software. The reason of this choice is the existence of a long legacy software, since a first, very sketchy, version was designed around 1995. mmalpha makes use of the polyhedral library polylib [START_REF]The Polylib Library[END_REF] and uses a parser and a pretty-printer written in C. The pip software is also used.

The interest of using mathematica lies in some of the functions available in this software, for example, the linear programming solver which is directly available.

B vhdl Translation of one Register Equation

The following vhdl code is the translation of a register equation of the form 𝑌 [𝑡, 𝑝] = CopyBMirrMovedIn[𝑡 -𝑝, 𝑝] that appears in the generation of the full Sherman-Morrison version, for 𝑁 = 13. It contains a Generate statement for the 𝑝 index, and another one for the 𝑖 index needed for an intermediate register named YReg. Inside the Generate statements, there is a sequential process needed to produce registers. The 𝑌 [𝑡, 𝑝] value is obtained through an assignment which is controlled by the value of the time counter.

14 2. 4 . 1

 1441 Onduleur triphasé 2 niveauxLe premier modèle à l'étude est un onduleur triphasé 2 niveaux[START_REF]The Polylib Library[END_REF]. Le circuit est représenté à la figure 2-1.

Figure 2 - 1 :Figure 1 .

 211 Figure 2-1: Schéma onduleur 2 niveaux En faisant la modélisation à l'aide de l'approche MNA, on observe que le modèle possède 11 noeuds et 5 sources de tension (V+, V-, Va, Vb, Vc). La matrice A sera donc une matrice de dimension 16×16 et le modèle sera le suivant :

14 11 noeuds et 5

 145 sources de tension (V+, V-, Va, Vb, Vc). La matrice A sera donc une matrice de dimension 16×16 et le modèle sera le suivant :

Figure 2 . 11 {Figure 3 .

 2113 Figure 2. The admittance matrix for the circuit of Fig. 1. The 𝑔 1 to 𝑔 6 values correspond to the admittance of the 6 switches 𝑆 1 to 𝑆 6 . Only these values are subject to a modification, when the status-open or closed-of the switches changes.

Figure 4 .

 4 Figure 4. Alpha code of the Sherman-Morrison algorithm

1 --Translation of the definition of Y 2 --Equation (register): 3 - 4 - 5 G81 18 (

 1234518 -Y[t,p] = CopyBMirrMovedIn[t-p,p] -Domain: {t,p | p+1<=t<=p+13; 1<=p<=13} counter -counterDelay >= p+1)

  . Each variable 𝑋 is scheduled by means of an affine-by-variable integer schedule of the form 𝑡 𝑋 (𝑧) = 𝛼 𝑋 (𝑧) + 𝛽 𝑋 . Causality is enforced by building an Integer Linear Program.

	Schedule may be multi-dimensional, but in our exam-
	ples, this possibility is not used.
	4. The alpha program is then transformed by applying
	to each variable a time-space change of basis. Each
	schedule is completed automatically by an affine allo-
	cation function in order to obtain a unimodular trans-
	formation (most often, this amounts to complete the
	schedule with a projection on 𝑛 -1 indexes of the vari-
	able). At the end of this step, all variables of the pro-
	gram are indexed by 𝑡 and 𝑝 (where 𝑝 represents a list
	of processor coordinates).
	5. The time-spaced alpha program is massaged into a
	multi-dimensional RTL description, by means of sim-
	ple transformations. After this step, equations are sep-
	arated into:
	• Combinational equations of the form
	X[t,p] = f(Y[t,p], ...)
	where 𝑓 is a function involving only combinational
	operators.
	• Simple connections, of the form X[t,p] = Y[t,p].
	• Register equations, X[t,p] = Y[f(t,p),g(t,p)].
	• Multiplexers, expressed using if and case expres-
	sions.
	• Control equations, involving boolean values.

When the time is multidimensional, which is not considered here, the notion of life-time is still applicable, with possibly several dimensions.

The vhdl code of the optimized version was "patched" to implement the direct access to the columns of the 𝐴 -1 matrix. In alpha, addressing an element of a variable through an index-say, writing X[A] where A is the result of an equation-is not direct, whereas the translation of such a statement in vhdl is obvious. The modification of the final code by hand concerned a few declarations and instructions. Obtaining the same code automatically is not, in theory, a difficulty, but it requires the mmalpha tool to be extended and checked.

Actually, without the calculation of the 𝜎 term, as it is supposed to be pre-calculated.

Acknowledgments

This work was funded by the Natural Sciences and Engineering Research Council of Canada grants, Prompt, Opal-RT, Hydro-Québec, CMC Microsystems, and the Research Chair in Signals and Intelligence of High-Performance Systems.

The authors would like to answer the referees for their insightful comments and suggestions to improve the clarity of this paper.

(counter -counterDelay <= p+13);

21

END GENERATE;


Pakedam Lare

Siyamak Sarabi

Claude Delpha

Demba Diallo

A PMaSynRM stator winding fault detection approach using an optimized PCA-based EWMA control scheme

Keywords: Fault detection, PMaSynRM Inter-turn shortcircuit fault, PCA, EWMA, optimization

This paper aims to present a statistical methodology for inter-turn short-circuit fault detection in a PMaSynRM for electrified vehicles in closed-loop operation. The proposed methodology is based on the combination of the Principal Component Analysis (PCA) and Exponentially Weighted Moving Average (EWMA) whose performance strongly depends on its control parameters. The False Alarm (PFA) probability is minimized to tune the control parameters correctly. Additionally, we develop a strategy based on comparing the eigenvalues of the data in the PCA subspaces to automatically select the statistical index allowing the best fault detection with the minimum missed detection probability (PMD). A hybrid Finite Element-Analytical model of a PMaSynRM has been simulated to generate synthetic data. For seven fault severities and three different loads, the performance is evaluated. We show that for all operating and fault conditions, fault detection is ensured in the worst case with P M D < 0.02 and P F A < 0.055. The proposed methodology outperforms traditional approaches, and its efficiency is proven.

I. INTRODUCTION

The reliability of electric motor drive systems in electric vehicles (EVs) has attracted considerable interest in recent years. Among the electric machines faults, inter-turn shortcircuits are the most common ones [START_REF] Gandhi | Recent advances in modeling and online detection of stator interturn faults in electrical motors[END_REF]. They are usually the cause of several other winding defects [START_REF] Nandi | Condition monitoring and fault diagnosis of electrical motors-a review[END_REF]. Then, for security reasons, early inter-turn short-circuit fault detection is very important. To generate the motor data to analyse for fault detection, a hybrid Finite Element-Analytical model of a three-phase Permanent Magnet-Assisted Synchronous Reluctance Motor (PMaSynRM) developed in [START_REF] Lare | Stator winding inter-turn short-circuit and air gap eccentricity fault detection of a permanent magnet-assisted synchronous reluctance motor in electrified vehicle[END_REF] and [START_REF] Lare | Modelling of a pma-synrm for the detection of inter-turn short-circuit[END_REF] has been developed. It has been simulated in closed-loop with a two-leve inverter, and validated under different operating conditions.

In recent inverter-fed motor applications, conventional motor monitoring tools face challenges due to high noise level, the dynamic changes that affect the excitation frequencies and closed-loop action that could affect the fault signatures [START_REF] Briz | High-frequency carrier-signal voltage selection for stator winding fault diagnosis in inverter-fed ac machines[END_REF]. To meet the need for comprehensive failure diagnosis tools, valid for all applications under different operating conditions, statistical analysis tools are an alternative to conventional methodologies [START_REF] Ayhan | Statistical analysis on a case study of load effect on psd technique for induction motor broken rotor bar fault detection[END_REF]. For multivariate systems, the Principal Component Analysis (PCA) has been widely used because of its good performance in highlighting variables' correlation and its no prior knowledge requirement [START_REF] He | Process fault detection and diagnosis based on principal component analysis[END_REF]. However, conventional PCA-based fault detection solutions using the fault indices Hotelling T 2 and Squared Predicted Error (SP E), are less suitable to detect mild changes [START_REF] Wang | Sensor validation and reconstruction for building central chilling systems based on principal component analysis[END_REF]. This limitation leads to the adoption of the Exponentially Weighted Moving Average (EWMA) monitoring chart where the state of the process at a point depends on the exponentially weighted average of all prior data. The sensitivity of the EWMA control procedure to gradual faults in the process highly depends on the choice of its tuning parameters [START_REF] Harrou | Amalgamation of anomalydetection indices for enhanced process monitoring[END_REF]: the weighting factor η and the control limits' parameter L.

In this paper, a methodology to optimally tune the parameters η and L is proposed. We show how to calibrate these parameters using an optimization scheme that minimizes the probability of false alarm (PFA). This methodology is illustrated in the case of inter-turn short-circuit fault detection for a PMaSynRM. Additionally,. The issue of noise influence in PCA feature space selection is addressed using an eigenvalue comparison test. Then, the most efficient framework and index can be automatically selected leading to the reduction of the probability of missed detection (PMD). The rest of the paper is organised as follows. In section II, after recalling the basics on PCA and EWMA for fault diagnosis, the proposed methodology is presented and the optimization procedure described. Section III presents the results and performance of the method. Section IV concludes the paper.

II. FAULT DETECTION METHODOLOGY

Let us consider N observations for n measured variables arranged in a data matrix X ∈ R N ×n , with zero mean and unit variance. In this matrix, each vector is written as

x T i = [x i,1 , . . . , x i,N
] where x i,j is the jth observations of the ith variable. In this approaches, the collected data in healthy operating conditions are used to build a reference model. Then the elements of this model are used to evaluate the faulty data.

The fault detection methodology proposed in this paper is based on the improvement of the limitations shown in conventional approaches using PCA and EWMA. These two approaches are first described and the improved optimized solution is further derived.

A. Principal Component Analysis for fault detection

Basically, PCA is used to reduce data dimensionality by grouping correlated variables in a set of new uncorrelated ones [START_REF] Jolliffe | Principal component analysis: a review and recent developments[END_REF]. Then, the data matrix X is subject to a linear transformation matrix T to express the variation in the observations such as T = XP, with T ∈ R N ×n the principal component matrix and P ∈ R n×n the eigenvectors matrix containing the eigenvalues λ i of the correlation matrix Σ of X such as:

Σ = P∧P T (1)
Where ∧ = diag(λ 1 , . . . , λ n ) is a diagonal matrix with diagonal elements in decreasing magnitude order. The reduced l-dimensional space of the data is obtained by retaining only the principal components that correspond to the l-highest eigenvalues of the covariance matrix. They represent the feature space where most of the pattern in the data are represented (generally up to a given percentage of cumulated variance).

The remaining components with the smallest eigenvalues represent the noise contribution. Once the number of components l to retain is determined, the eigenvectors matrix P and the principal component matrix T are partitioned into the form:

P = Pl P n-l T = Tl T n-l (2)
X is the principal part of the data explained by the l first eigenvectors and the residual part X is explained by the residual components [START_REF] Tharrault | Fault detection and isolation with robust principal component analysis[END_REF]:

X = Pl PT l X X = P n-l P T n-l X (3)
They respectively lead to the principal and residual subspaces used for fault diagnosis.

A process monitoring with PCA uses Hotelling's T 2 and Q statistics also called Squared prediction error (SP E) to detect abnormal behaviours. The Hotelling's T 2 and SP E represent the variability in the Principal Component Subspace (PCS) and the Residual Subspace (RS), respectively. T 2 can be expressed by using the estimated value Σ of the correlation matrix as follows [START_REF] Engelen | A comparison of three procedures for robust PCA in high dimensions[END_REF]:

T 2 = XT Σ-1 X (4)
Under normal conditions if the data set is multivariate Gaussian distributed, the T 2 can be approximated by a χ2 distribution with l degrees of freedom and a significance level (1-α).

The system operation is considered healthy if:

T 2 ≤ χ 2 l,1-α (5)
A variation of the variables' correlation indicates an abnormal behaviour. Under this condition, the sample matrix X increases its projection to the RS and the magnitude of X reaches abnormal values compared to those obtained during healthy conditions. SP E is the magnitude of X and is written as:

SP E = X 2 (6)
The process is considered healthy if the SP E statistic is under its control limit which is expressed as follows:

SP E ≤ (μ + σz 1-α ) 3 (7)
with z 1-α as the (1α) significance level of a Gaussian distribution Φ. μ and σ are respectively the estimated mean and standard deviation of SP E 2/3 [START_REF] Engelen | A comparison of three procedures for robust PCA in high dimensions[END_REF].

Unfortunately the process monitoring using T 2 and SP E is not always effective as they only use the current observation information for decision making. Hence, these detection indices are relatively insensitive to minor changes in the process variables, and thus may result in missed detections. To overcome these limitations, an alternative approach has been developed, where PCA is used as a modelling framework for fault detection using an EWMA control scheme [START_REF] Harrou | A statistical fault detection strategy using PCA based ewma control schemes[END_REF].

B. PCA-based EWMA control scheme for fault detection

In the EWMA control scheme, the moving average is calculated by multiplying the historical observations with weight coefficient that exponentially decays with time. The EWMA decision statistic criterion is then described as follows:

z t = η × x i,t + (1 -η)z t-1 (8)
where η is a weight parameter, with 0 < η ≤ 1 , and x i,t is the value of the monitored variable i at time t. The starting value x i,0 is set equal to the mean µ 0 of the process in-control.

Generally, smaller values of η increase the chart's sensitivity to smaller shifts in the process mean, while larger values of η increase its sensitivity to larger shifts. The standard deviation of z t is approximated as follows:

σ zt = η 2 -η (9)
The EWMA control scheme declares an anomaly when the value of z t falls outside of the interval between the upper and lower bounds of the control limits . The upper and lower control limits, denoted as U CL and LCL are set as:

U CL = µ 0 + Lσ zt LCL = µ 0 -Lσ zt ( 10
)
Where L is a multiplier of EWMA standard deviation σ zt . As stated in Section I, the sensitivity of the control procedure to small or gradual faults highly depends on the settings of L and η. The fault detection using PCA-based EWMA scheme involves the evaluation of these residuals from the PCA framework using a statistical EWMA to set the control limits U CL and LCL. It allows to determine if the system is healthy or faulty. If the EWMA statistic exceeds any of the control limits, this is an indicator of an abnormal (faulty) behavior [START_REF] Harrou | A statistical fault detection strategy using PCA based ewma control schemes[END_REF].

C. Optimized proposed fault detection scheme

The proposed solution for this study has to tackle two main limitations of the EWMA control scheme: (1) the selection of the correct PCA feature space, (2) the tuning of the weight and the control limit parameter.

1) In the conventional PCA-based EWMA fault detection algorithm, one major limitation using the EWMA monitoring scheme is in the exclusive use of the residual components of the PCA model. Indeed, this is valid in ideal systems where the presence of faults affects mostly the components in the residual subspace, but this is not always true. In some cases, the system modelling errors could be projected onto the residual subspace making it less sensitive to the fault and the principal subspace more sensitive. So, it is necessary to analyse the information in both subspace as a whole set to obtain more accuracy in the fault diagnosis results.

In our work, we propose to apply the EWMA monitoring scheme to either the principal and residual subspace depending on the evaluation of the eigenvalues of the test data in the PCA subspaces. We compare the smallest eigenvalue in the principal subspace to the highest one in the residual subspace. The higher the ratio between these aforementioned eigenvalues is, the more the fault signature is found highly represented in the residual component and vice versa. The correct subspace can then be automatically selected. 2) Usually, the EWMA control limit parameters L and η are chosen depending on the type of system under monitoring.

In most of applications, the EWMA standard deviation amplifier L is set equal to 3 while the weighted parameter η is set depending on the importance one wants to give to the current observation in the calculation of the EWMA.

To accurately tune these parameters, we propose to set them from the solution of an optimization problem. For that purpose we show how to obtain them following the minimization of the probability of false alarm (P F A) on healthy data. This probability is a useful criterion to set a preliminary detection threshold without any a priori knowledge on the faulty data. All the data that are above this threshold will be falsely considered as faulty (false positive). The others are truely considered as healthy (true negative). Then the P F A is given as :

P F A = N F P N F P + N T N (11)
Where N F P is the number of false positive values and N T N is the number of true negative ones. The problem in our case is to find the couple of (η, L) that minimizes the P F A. The fact that the conventional P F A parameters N T N and N F P are only defined by crossing the two EWMA control limits makes the problem divergent. To get a convergent problem we define a new objective function relating the distance between the EWMA observations z i and their control limits U CL and LCL as follows:

min η,L (z max (η j , L j ) -U CL(η j , L j )) 2 min η,L (z min (η j , L j ) -LCL(η j , L j )) 2 s.t. 0 < η j ≤ 1, 0 < L j ≤ 3 (12)
Where (η j , L j ) are the potential solutions, z max , z min are respectively the maximum and the minimum of the EWMA observations. The multi-objective function can be reformulated into a mono-objective one such as:

min η,L F obj (η j , L j ) s.t. 0 < η j ≤ 1, 0 < L j ≤ 3 (13)
Where:

F obj (ηj, Lj) = (zmax(ηj, Lj) -U CL(ηj, Lj)) 2 + (zmin(ηj, Lj) -LCL(ηj, Lj)) 2 (14)
The proposal is described in the flowchart displayed in Figure 1.

III. RESULTS AND DISCUSSIONS

In this section, the PCA-based EWMA methodology is used to detect stator inter-turn short-circuit fault in a PMaSynRM.

A. Data generation

To assess the performance of the proposed fault detection scheme, the model of the PMaSynRM presented in [START_REF] Lare | Stator winding inter-turn short-circuit and air gap eccentricity fault detection of a permanent magnet-assisted synchronous reluctance motor in electrified vehicle[END_REF] and [START_REF] Lare | Modelling of a pma-synrm for the detection of inter-turn short-circuit[END_REF] is used. The motor characteristics and simulation parameters are displayed in Table I. The PMaSynRM winding under an inter-turn short-circuit in phase a in shown in Figure 2.

Under this fault condition, the phase a current i a is divided into i s the current in the defective part of the winding and i f the current in the inter-turn fault's contact branch. L a1 and L a2 represent the inductances of the healthy and faulty sub-coils of the phase a winding respectively. In one hand, M a1b and M a1c denote the mutual inductances between the healthy subcoil L a1 and the coils L b , and L c respectively. On the other hand M a1a2 , M a2b and M a2b denote the mutual inductances between the sub-coil L a2 and the healthy sub-coil L a1 , and coils L b , and L c respectively.

Three different mechanical load conditions (Γ 1 = 18N m, Γ 2 = 134N m, Γ 3 = 260N m) have been applied to the machine during 3s each, with a controlled rotating speed at 1500rpm. The generated data include the motor stator phase currents i a , i b and i c . The faults which are changes in the values of the stator resistances and inductances are introduced after 1s of healthy state for each load case. The duration of the fault is also 1s. To assess the abilities of the various fault detection of the methodology, seven different severity levels are considered in this study. The fault severity levels are ranked

Data collection

Healthy and faulty data

Training data

Healthy data

Test data

Faulty data

Normalisation of training data

Normalisation of test data

PCA model Projection on the PCA spaces

Compute T² and SPE

Compute EWMA observations z i (η j ,L j ), z max (η j ,L j ), z min (η j ,L j ), UCL(η j ,L j ) and LCL(η j ,L j )

for T² and SPE

Set of initial population

(η j ,L j ) from to 2 to 8. This rank is related to the number of turns shortcircuited in the phase a of the winding. As an example, the time series stator's current for the phase a of the PMaSynRM drive is shown in Figure 3. The fault occurrence with severity level 2 and load change (Γ 1 , Γ 2 , Γ 3 ) are displayed highlighting the healthy and faulty current information.

Compute

In our case studies, using the time series currents of the three phases, the statistical moments are first calculated for consecutive independent windows of the current signal samples containing at least one signal period. The goal of these features is to highlight the quantitative insights of these currents and be used as the variables for the fault diagnosis study. Among these

Fault introduced

Fault ended

Fault ended

Fault ended

Fig. 3: Time series stator current in phase a for the healthy and faulty motor (level 2) for the three loads features, the variance appears to be the best representative one to characterize the considered fault signature. Then the variances, σ 2 a , σ 2 b , σ 2 c of the phase currents in healthy and faulty conditions have been calculated and organised in the data matrix to be used as the input of our fault detection procedure. The following results highlight the performances of our proposal.

B. Simulation results

As stated in Figure 1, the procedure is decomposed in two parts: one using the healthy data for modelling the process and the other considering the faulty data for the fault detection. With the healthy data, the process is trained to construct the PCA model. In this study, the threshold of cumulative variance is settled to 75%. In our PCA model, this results in retaining the first two components for the principal subspace and the third one for the residual subspace. This principal subspace captures in the first two principal components 43.31% and 32.94% for Γ 1 , 75.62% and 14.36% for Γ 2 , 53.82% and 24.91% for Γ 3 (see Figure 4).

With the faulty data, the fault detection ability of the proposed PCA-based EWMA fault detection algorithm is assessed. All the seven severity levels and the three load levels are considered. The conventional PCA-based EWMA fault detection algorithm with the application of T 2 and SP E In this Figure 5 the dashed lines represent U CL and LCL. The results using PCA-based EWMA on the SP E statistic, show that it failed to detect some faults and produce a high number of false alarms. The results using PCA-based EWMA on the T 2 statistic allow to detect the fault with few false alarms. This result is because the conventional PCA-based EWMA with the arbitrary choice of the parameters η and L is not always sensitive to the fault due to the control limits that are not particularly calibrated with the considered healthy data.

Figure 6 shows the overall performance in terms of probability of missed fault detection P M D and probability of false alarm P F A for our case studies of three load levels and all fault severity levels. Usually, for fault detection, the acceptable P M D and P F A are typically set to 0.02 and 0.05 respectively [START_REF] Weber | Increasing effectiveness of model-based fault diagnosis: A dynamic bayesian network design for decision making[END_REF]. In Figure 6, these typical values are highlighted with dashed lines. One can notice that the conventional PCA-based EWMA fails to provide acceptable P F A using T 2 and SP E. The P M D performance for low severity levels are good only when using T 2 for the highest load conditions. These results are also acceptable when using SP E for the lowest load conditions. Indeed, the fault signatures are not always high enough in the residual subspace to allow the detection in the presence of noise.

To cope with these limitations, our proposed methodology include the optimization tuning of η and L for a minimized P F A value through an objective function. Figure 7 shows the surface of the objective function using T 2 and SP E for a The overall performance with the optimized couples (η, L) defined for each load condition is shown in Figure 8. Compared to the conventional approach, the performance is improved but show a sensitivity to the load value and fault severity level depending on the chosen detection index. The last step of the proposal concerns the automatic selection of the detection index by using the ratio of selected eigenvalues in the PCA framework. In our case study the reference of the ratio is set equal to 2 to set the most efficient fault detection index. A ratio less than 2 sets the T 2 as detection index while a ratio greater than 2 sets the SP E as the detection index.

Figure 9 shows the results of the fault detection of the proposed optimized PCA-based EWMA scheme with the automatic selection of the detection index in the case of fault severity level 2, and for the three load conditions. The overall performance of the proposed methodology is shown in Figure 10. This result highlights that the proposed methodology can reach good fault detection performance as the P M D and P F A almost comply with the usual typical values. The only exception is for fault severity 8 and load 2. In that case, the P F A is slightly higher (P F A = 0.052) than the typical value without being critical. Most of the False alarms in our case studies are due to the transient between the faulty and healthy conditions. In this paper, a fault detection approach is proposed for inter-turn short-circuit monitoring in the electrical drive. A hybrid Finite Element-Analytical model of a PMaSynRM has been simulated to generate the stator phase currents for analysis. Three load conditions and seven fault severities are considered. The proposal is based on PCA and EWMA that are efficient techniques for fault diagnosis. Nevertheless, these techniques suffer from two major limitations: (1) the tuning operation of the control parameters to reach the best performance, (2) the selection of the most relevant transformation subspace with the adequate detection criteria.

This paper proposes solutions to cope with these limitations and an optimized procedure with the automatic selection of the most efficient detection criteria. The paper presents how to optimize the EWMA weighting factor η and control limit parameter L with the constraint of minimizing the probability of false alarm. Additionally, we show that the comparison of the eigenvalues of the test data in the PCA subspace allows to automatically choose the most efficient subspace framework and adequate fault index. The simulation results show encouraging performance. Compared to the traditional approaches, the advantage of the proposed optimized PCA-based EWMA fault detection method is demonstrated whatever the load condition and for all the inter-turn short-circuit severity levels.

In future works, the application to real data will be studied, and the influence of disturbances will be included. The reliability of the decision will also be evaluated in terms of mean delay for detection and mean time between false alarms. Moreover, the performance of the proposed approach to detect other faults in the electrical drive will be investigated.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Optimized PCA-based EWMA methodology

Fig. 4 :

 4 Fig. 4: Variance captured by each principal component

Fig. 5 :

 5 Fig. 5: Conventional PCA-based EWMA fault detection

6 PFA for T 2 F a u lt le v e l 2 F a u lt le v e l 3 F a u lt le v e l 4 F a u lt le v e l 5 F a u lt le v e l 6 F a u lt le v e l 7 F a u lt le v e l 8 0F a u lt le v e l 2 F a u lt le v e l 3 F a u lt le v e l 4 F a u lt le v e l 5 F a u lt le v e l 6 F a u lt le v e l 7 F a u lt le v e l 8 0Fig. 6 :Fig. 7 :

 622345678234567867 Fig. 6: Performance of the conventional PCA-based EWMA fault detection

F a u lt le v e l 2 F a u lt le v e l 3 F a u lt le v e l 4 F a u lt le v e l 5 F a u lt le v e l 6 F a u lt le v e l 7 F a u lt le v e l 8 0F a u lt le v e l 2 F a u lt le v e l 3 F a u lt le v e l 4 F a u lt le v e l 5 F a u lt le v e l 6 F a u lt le v e l 7 F a u lt le v e l 8 0Fig. 8 :

 234567823456788 Fig. 8: Performance of the improved PCA-based EWMA with two fault indices

Fig. 9 :

 9 Fig. 9: Fault detection with the improved PCA-based EWMA with automatic detection index selection

F a u lt le v e l 2 F a u lt le v e l 3 F a u lt le v e l 4 F a u lt le v e l 5 F a u lt le v e l 6 F a u lt le v e l 7 F a u lt le v e l 8 0F a u lt le v e l 2 F a u lt le v e l 3 F a u lt le v e l 4 F a u lt le v e l 5 F a u lt le v e l 6 F a u lt le v e l 7 F a u lt le v e l 8 0Fig. 10 :

 2345678234567810 Fig. 10: Performance of the improved PCA-based EWMA automatic detection index selection

UCL and LCL or SPE

	Compute and Minimize F
	(η j ,L j ) optimized	(η j ,L j ) optimized
	for T²	for SPE
	Compute UCL
	and LCL or T²

obj (η j ,L j ) for

T² and SPE Evaluate the eigenvalues ratio Determine fault if z i >UCL or z i <LCL Determine fault if z i >UCL or z i <LCL

	No	Yes
	Compute T²	Compute SPE
	Compute EWMA	Compute EWMA
	observations z i	observations z i

TABLE I :

 I The PMaSynRM characteristics

	Motor characteristics
	Number of poles	8
	Number of PMs per pole	5
	Rated power [kW ]	208
	Maximum operating speed [rpm]	14000
	Number of turns	8
	Stator phase winding resistance value at 20 • [mΩ]	6.1
	Simulation parameters
	Permanent Magnet flux amplitude [W b] 6.24 ×10 -2
	Moment of inertia [kgm 2 ]	0.0357
	Input DC voltage amplitude [V ]	412

Authorized licensed use limited to: Centrale Supelec. Downloaded on January 17,2024 at 10:39:13 UTC from IEEE Xplore. Restrictions apply.