Erwan Pierre
email: erwan.pierre@elmy-energie.fr.

Lorenz Schneider
email: schneider@em-lyon.com.

Arthur Boidevézi

Rüdiger Kiesel

Jules Lafritte

Cassio

Intermittently coupled electricity markets

Keywords: Electricity Markets, Interconnectors, Market Coupling, Spread Options, Regime Switching JEL: C51, C63, D44, G13, Q41

. We describe the model estimation procedure in detail, and compare model and market prices of European spread options. We observe a dramatic paradigm shift in our data set at the end of the summer of 2021, and show that this shift has a strong effect on the model parameters. We also see that the reliable pricing and trading of spread options becomes problematic in such a volatile and uncertain market environment.

Introduction

European interconnected electricity markets

From the very beginning of its efforts to liberalise the European electricity markets in 1996, the European Commission has focused on ensuring the development of a secure and integrated market. In line with this objective, it set up an Expert Group to provide specific technical advice. The European Commission Expert Group (2017) states in its first report that "Electricity interconnectors are the physical component of making this market truly European by connecting Member States' networks offering capacity for electricity trade, improved security of supply and allowing integration of the rapidly-growing share of renewable electricity production." The European Commission set the objective of increasing the electricity interconnection target to 15% by 2030, this ratio being defined as import capacity over installed generation capacity in a given Member State.

Since then, the European Commission Expert Group (2018) has made significant changes to this ratio to take into account the increasing share of renewables in the European energy mix. And in its Market coupling and the role played by interconnectors have also been the subject of several studies in recent years. [START_REF] Parisio | Electricity prices and cross-border trade: Volume and strategy effects[END_REF] introduce a setup with equilibrium bid functions in individual domestic markets. They study the way these functions are affected when cross-border trading is managed using the implicit auction method, and analyse the effect of this mechanism on price convergence. With a particular focus on Italy, [START_REF] Pellini | Measuring the impact of market coupling on the Italian electricity market[END_REF] evaluates the impact on the Italian electricity market of replacing an explicit auction mechanism with an implicit one, i.e., with market coupling, and shows that the effects on welfare and social surplus are likely to be positive. [START_REF] Cartea | How much should we pay for interconnecting electricity markets? A real options approach[END_REF] take a real options approach to determine the value of an interconnector. They study five pairs of neighbouring European markets, determine the values of the interconnectors for each pair over a one-year period in their setup, and find that jumps in price spreads can have a significant impact on these values. [START_REF] Füss | Electricity spot and derivatives pricing under market coupling[END_REF] employ a fundamental multi-market model to study electricity spot and derivatives pricing under market coupling, and analyse how the behaviour of these prices depends on different cross-border allocation schemes. [START_REF] Kiesel | Structural models for coupled electricity markets[END_REF] calibrate their structural model to the French-German market area and study the effect of market coupling on electricity prices and the value of power plants. They find that market coupling can lead to lower futures prices in all affected markets, even the cheapest one, and can also result in volatility spillover effects. [START_REF] Pircalabu | A regime-switching copula approach to modeling day-ahead prices in coupled electricity markets[END_REF] combine regimeswitching with a copula model for a pair of electricity prices. Each hour is considered independently of the others, and pricing is performed using a Monte Carlo simulation. Their empirical study is based on four pairs of market areas, for which they find significant evidence of tail dependence in the interconnected areas considered. [START_REF] Parisio | Market coupling between electricity markets: Theory and empirical evidence for the Italian-Slovenian interconnection[END_REF] focus on the market coupling between Italy and Slovenia, and assess the degree of integration between these markets arising because of this policy. They find that despite this existing link, the two markets are still far from being strongly integrated. [START_REF] Cartea | Speculative trading of electricity contracts in interconnected locations[END_REF] analyse an investor's optimal strategy when trading electricity contracts in two interconnected locations. They show that as his aversion to ambiguity increases, his trading activity, and consequently his inventory exposure, decreases in both locations. [START_REF] Sønderby | Modelling the joint behaviour of electricity prices in interconnected markets[END_REF] introduce regime-switching in their model for potentially coupled electricity prices using a latent uni-variate process that allows for seasonal behaviour. The parameters of this process are estimated from given market data via particle filtering. They show that their model can be used for the analytical pricing of forwards and spread options, and illustrate its performance in a case study of the French and German markets. [START_REF] Simon | A stochastic study of carbon emission reduction from electrification and interconnecting cable utilisation[END_REF] consider the new NordLink cable between Germany and Norway, and analyse the effect that the exchange of power production between two countries can have on carbon emissions. They conclude that the cable can be used as an effective way to reduce carbon emissions by exporting surplus German renewable production to Norway on windy days, and importing surplus Norwegian energy production to Germany on non-windy days.

In this article, we take the [START_REF] Cartea | How much should we pay for interconnecting electricity markets? A real options approach[END_REF] modelling framework as our point of departure, and set out from there to take the aforementioned developments in market coupling and interconnectors into account. One such development is the finer granularity in electricity markets. [START_REF] Cartea | How much should we pay for interconnecting electricity markets? A real options approach[END_REF] consider the means of day-ahead prices (also called spot prices) over peak and off-peak hours when pricing interconnector options. In practice, however, these options are strips of individual options for each hour of the day. This distinction is important, because the dynamics of the hourly spreads are quite different from the dynamics of the spreads of mean spot prices. In fact, due to the physical constraints of the connection, the hourly spread is either equal to zero if there is enough capacity to allow the same price in the two locations, or different from zero if the lines are overloaded. To account for these network constraints, we propose a formulation of the problem with regime-switching for the underlying hourly spreads. We note that our reduced-form model is highly tractable and allows for robust parameter estimation. As we will show, prices of interconnector options can be computed in closed-form. Another even more recent development is the dramatic paradigm shift in European electricity markets that began in the summer of 2021. Drawing on our data set, which starts in 2017, we illustrate the strong effect of this shift on the model parameters. Finally, we consider the issues of reliably pricing and trading spread options in such a volatile and uncertain market environment.

Main contributions

Our main contributions are:

• We increase the resolution from daily (peak, off-peak) electricity prices to hourly prices, and add regime-switching to reflect coupled and decoupled hours.

• We perform a simulation exercise to demonstrate that our regime-switching model successfully captures the important stylised features of intermittently connected markets.

• We show that our model can adapt to the new market paradigm of high prices, uncertainty and volatility.

• We find that the pricing of spread options works well in a "normal" market environment. However, in a highly volatile environment, both auction and model ex ante spread option prices may turn out to lie far from ex post realised payoffs.

Structure of the article

The remainder of the article is organised as follows. In Section 2, we describe the current European interconnector market. In Section 3, we present our model, which is based on the Cartea and González-Pedraz (2012) model and a 2-state Markov process that describes regime-switching between coupled and decoupled periods in the market. Section 4 explains how to estimate the model parameters for given hourly power prices from various European electricity markets. In Section 5, we show that our model can adapt to the new market paradigm that began in early autumn, 2021. Section 6 concludes the article. In Appendix A, we review Fourier transform methods for the pricing of European spread options.

The European Interconnector Market

An Interconnected European Electricity Market

The interconnection market in Europe allows market participants to buy or sell flows of energy between two locations in order to exploit price differentials between locations. These locations are delimited in Europe by the different Transmission System Operators (TSOs) and can be a country or a larger bidding area. For example, France (FR) has interconnections with Belgium (BE), Switzerland (CH), Germany (DE), Spain (ES) and Italy (IT).

Before the introduction of this market, a trading member had to first reserve cross-border capacity, before transporting the purchased electricity in a second step. The interconnection market simplifies the process via implicit auctions in which market participants simply bid for electricity on the exchange. This auction market comprises both intra-day products and long-term options. When a participant buys an option in the DE→FR interconnection, for example, he buys the right to transit energy for each hour of the period considered from Germany to France. When the flux is financial, the payoff of such an option is:

Payoff = h∈H max S A (h) -S B (h) -K, 0 , (1)
where K is the strike price, S c (h) is the spot price for each country c = A, B at a given hour h, and H is the set of all hours in the period. The price of such an option is a function of the dynamics of the underlying spread. Figure 1 displays the hourly spread S F R,DE = S F R -S DE used to value the DE→FR interconnection. We observe that during many hours in this period, electricity prices were higher in France than in Germany, leading to positive price spreads.

When electricity is sent from one country to another, the ensuing spread can be denoted in two ways that must not be confused: DE → FR means that electricity is sent from Germany (DE) to France (FR), which means that electricity is bought in DE and sold in FR, which leads to a spread of S F R,DE = S F R -S DE . In the other direction, FR → DE means that electricity is sent from France to Germany, which means that electricity is bought in FR and sold in DE, leading to the opposite spread of S DE,F R = S DE -S F R .

Market Coupling and Decoupling

In Europe, the day-ahead price is fixed by an algorithm called Euphemia and powered by European electric power exchange EPEX Spot. The algorithm allows us to take into account the Price Coupling Regions (PCR) by optimising the maximal welfare with the interconnection constraints. The efficiency of market coupling is proven by an increasing price convergence between market areas. As a consequence, when the interconnection lines between two countries for a given hour are not saturated, the price is the same for both, which means that the spread is zero. In this case we say that the market is coupled. In contrast, when there is insufficient capacity and prices differ, we say that the market is decoupled.

If we take Figure 1 and zoom in over a shorter period, for example May 2018, as shown in Figure 2, we observe that there are many occasions when the hourly spread is zero. In this example, 20% of the DE and FR hours are coupled, with an ensuing price spread of zero.

J o u r n a l P r e -p r o o f

Journal Pre-proof market price spreads with neighbouring countries for 24 May 2022 from EPEX Spot. We can see that at the beginning of the day, in hours 1 to 3, all German spreads are different from zero, meaning that Germany is completely decoupled from its neighbours. Then, for most of the day, some countries are coupled with Germany, while others are not. Finally, in hours 21 to 23, all German spreads are equal to zero, meaning that Germany is completely coupled with its neighbours.

The situation can change if we consider daily rather than hourly price spreads. To illustrate the different behaviour of hourly and daily spreads in the market, Figure 3 shows the daily spread in the first half of 2018. Due to the averaging of the hourly prices over the day, and since there are virtually no days on which all 24 hours are coupled, the coupled regime disappears, and only the decoupled regime is observed. The periods over which price spreads are calculated when determining the payoff of an interconnection option are therefore an important contractual specification.

To understand the importance of these periods, let us take a closer look at the price spreads observed in the market between Germany and Poland shown in Table 1. The sum of the spreads over the eight peak hours from 9 to 20 in column DE→PL is -121.56 (in EUR/MWh), or -10.13 per hour, on average. For the remaining off-peak hours from 1 to 8 and 21 to 24, the sum of the spreads is -155.74, or -12.98 per hour. Finally, over the whole day, the total base spread is -277.30, or -11.55 per hour. In all three cases, i.e., peak, off-peak and base load, the payoff of an interconnection option with strike K = 0 is equal to zero, since the three spreads are negative. This changes, however, if we look at an interconnection option whose payoff is based on each hour individually. The spreads for the hours 1 to 3 and then again 14 to 20 are positive, and when we calculate the payoff of this hourly option using equation ( 1), we obtain a positive payoff of EUR 186.02. Table 2 compares these different option payoffs for the couple DE→PL just discussed, as well as the other country couples of Table 1. The couple DE→PL stands out, of course, because of its pronounced alternation between positive and negative price spreads.

Electricity Spreads and Spread Options

Let S A (t) and S B (t) at time t be electricity spot prices in locations A and B, respectively, and define their spread as:

S A,B (t) = S A (t) -S B (t)
namely, the difference between these two prices. Let r be the continuously-compounded risk-free interest rate. The model price of a single call option with strike K A,B at an hour T is given by the discounted expected payoff:

C A,B (t, T, K A,B ) = e -r(T -t) E t max(S A,B (T ) -K A,B , 0) ,
and the value of the monthly option is then given by:

C A,B M (t, K A,B ) = TN T =T1 C A,B (t, T, K A,B ),
where T 1 ,..., T N are the hours of delivery in month M . In practice, the strike K A,B is usually set to zero. However, we consider the more general case of any strike value here. Non-zero values of the strike could be taken to account for possible interconnection costs, for example.

As usual in electricity markets, option prices are quoted per MWh. To obtain the actual price of the option, this price needs to be multiplied by the number of hours in the month.

If the characteristic function of the spread S A,B (T ) is available, but the probability density is not, option valuation theory allows us to use the Fourier-transform (see [START_REF] Carr | Option valuation using the Fast Fourier Transform[END_REF]) to numerically compute the price of such a spread option.

3 An Interconnector Model with Regime Switching

The Cartea-González-Pedraz Model

We begin the description of the interconnector model by recalling the main features of the Cartea and González-Pedraz ( 2012) model. They propose the following model for electricity spreads under the physical probability measure P. Let (Ω, A, P, F) be a filtered probability space. Let S A (t) and S B (t) be electricity spot prices at time t in locations A and B, respectively, and define their spread as:

S A,B (t) = S A (t) -S B (t), (2)
i.e., the difference between these two prices. [START_REF] Cartea | How much should we pay for interconnecting electricity markets? A real options approach[END_REF] assume that this spread follows an arithmetic model given by:

S A,B (T ) = f (T ) + X(T ) + Y (T ). ( 3
)
The seasonality function f describes a deterministic seasonal pattern. The process X is an Ornstein-Uhlenbeck process given by:

X(T ) = X(t)e -α(T -t) + T t e -α(T -s) σ(s)dB(s) (4)
with mean-reversion speed α > 0 and volatility σ > 0, and where B = (B t ) t≥0 is a standard F-adapted Brownian motion under P. The process Y is a zero-mean-reverting pure jump process given by:

Y (T ) = Y (t)e -β(T -t) + T t e -β(T -s) dJ + (s) + T t e -β(T -s) dJ -(s). (5)
with mean-reversion speed β > 0 and jump-intensities λ + , λ -≥ 0. In the general formulation of the model, these jump-intensities can be time-dependent. Here, for simplicity and better tractability, we assume that these parameters are constant.

The compound Poisson processes J + , J -are given by:

J + (t) = N + (t) n=1 j + n , (6)
J -(t) = N -(t) n=1 j - n , (7)
where N + (t), N -(t) are Poisson processes with intensities λ + andλ -, respectively. The variables j + and j -denote positive and negative jumps, which are exponentially i.i.d. with parameters η + , η -> 0. The expected jump sizes are therefore 1/η + and 1/η -.

As is shown in [START_REF] Cartea | How much should we pay for interconnecting electricity markets? A real options approach[END_REF], the characteristic function ϕ T = ϕ S A,B (T ) of S A,B (T ) is given for u ∈ R by:

ϕ T (u) = exp iuh(T ) - 1 4α u 2 σ 2 (1 -e -2αT ) • η + -iue -βT η + -iu λ + /β • η -+ iue -βT η -+ iu λ -/β , ( 8
) with h(T ) = f (T ) + e -α(T -t) X(t) + e -β(T -t) Y (t).
Note that ϕ T has singularities at u = -iη + and u = iη -.

J o u r n a l P r e -p r o o f

Journal Pre-proof

Regime Switching via a Two-State Markov Process

Let S = {0, ..., m -1} be a finite state space for a Markov process M = (M t ) t≥0 in continuous time on a probability space (Ω, A, P). Let the m × m-matrix Q = (q ij ) , q ij ∈ R ∀i, j ∈ S, satisfy the following conditions:

1)

q ij ≥ 0 ∀i, j ∈ S with i ̸ = j; q ii ≤ 0 ∀i ∈ S, 2) j∈S q ij = 0 ∀i ∈ S, (rows sum to zero) 3) 0 < max i∈S |q ii | < ∞.
The matrix Q is called the infinitesimal generator of the process M .

The transition probability matrices P (t) are given for t ≥ 0 by the matrix exponential:

P (t) = e tQ = ∞ n=0 (tQ) n n! . ( 9
)
The stationary distribution is defined by πQ = 0. It also satisfies πP (t) = π for all t > 0.

Assume that m = 2, S = {0, 1}, i.e., that the Markov process M has the two states 0 and 1. In this case, we simplify the notation and refer to the off-diagonal elements of Q as q 0 and q 1 , so that q 00 = -q 0 , q 01 = q 0 , q 10 = q 1 , q 11 = -q 1 . We assume that 0 < q 0 , q 1 < ∞. The generator Q is given by: Q = q 00 q 01 q 10 q 11 = -q 0 q 0 q 1 -q 1 . Using Kolmogorov's backward equation [START_REF] Stroock | An Introduction to Markov Processes[END_REF]:

P ′ (t) = QP (t),
with:

P (t) = p 00 (t) p 01 (t) p 10 (t) p 11 (t)
it can be shown that:

p 00 (t) = q 1 q 0 + q 1 + q 0 q 0 + q 1 e -(q0+q1)t , (10)

p 01 (t) = q 0 q 0 + q 1 - q 0 q 0 + q 1 e -(q0+q1)t , (11)
p 10 (t) = q 1 q 0 + q 1 - q 1 q 0 + q 1 e -(q0+q1)t , (12)
p 11 (t) = q 0 q 0 + q 1 + q 1 q 0 + q 1 e -(q0+q1)t .

(13)

Note that we have p 00 (t) + p 01 (t) = 1 and p 10 (t) + p 11 (t) = 1, which is natural from a Markov chain point of view. Note also that for the sum of the diagonal probabilities we have the condition:

p 00 (t) + p 11 (t) > 1, (14)
and for the sum of the off-diagonal probabilities we have the condition:

p 01 (t) + p 10 (t) < 1. ( 15
)
Conversely, if the transition matrix P (t) is given for some t > 0 and satisfies the restrictions ( 14) and ( 15), then the generator Q can be calculated from it as follows. This will be useful for estimating the parameters of the Markov process M .

Adding equations ( 11) and ( 12) gives:

p 01 (t) + p 10 (t) = 1 -e -(q0+q1)t ,
and therefore:

q 0 + q 1 = - 1 t ln(1 -p 01 (t) -p 10 (t)). ( 16
)
Now, substituting ( 16) in ( 11) and ( 12), respectively, gives:

q 0 = - 1 t ln(1 -p 01 (t) -p 10 (t)) • p 01 (t) p 01 (t) + p 10 (t) ( 17
)
q 1 = - 1 t ln(1 -p 01 (t) -p 10 (t)) • p 10 (t) p 01 (t) + p 10 (t) (18)
The stationary distribution, satisfying πQ = 0, is given by π = ( q1 q0+q1 , q0 q0+q1 ).

An Interconnector Model with Regime Switching

We now combine the Cartea and González-Pedraz (2012) model of Section 3.1 with the Markov regimeswitching process of Section 3.2.

Let S 0 = (S 0,t ) t≥0 and S 1 = (S 1,t ) t≥0 be two stochastic processes defined over the same probability space as the Markov process M . Further, we assume that S 0 and S 1 are both independent of the process M . We define the regime-switching process S as:

S(t) = S 0 (t)I {Mt=0} + S 1 (t)I {Mt=1} . ( 19
)
The characteristic function ϕ ST = ϕ T of S(T ) is then given for u ∈ R by: = P(M T = 0)ϕ 0,T (u) + P(M T = 1)ϕ 1,T (u), where ϕ 0,T and ϕ 1,T are the characteristic functions of S 0 (T ) and S 1 (T ), respectively.

The probability vector π = (π 0 , π 1 ), with 0 < π 0 , π 1 < 1 and π 0 +π 1 = 1, of the stationary distribution of the Markov chain satisfies the relationship:

πP (t) = π, or, equivalently, πQ = 0.
For the generator Q as given above, we obtain:

π = q 1 q 0 + q 1 , q 0 q 0 + q 1 .
In particular, note that π does not depend on time t.

We now come to the specification of the two regimes in our model. As we described in Section 2, the market switches between two regimes in which the spread S A,B (t) = S A (t) -S B (t) between electricity prices in two locations A and B at time t, as defined in equation ( 2), is either equal to 0 or different from 0. These are our regimes 0 and 1. We therefore specify that the process S 0 is constant and equal to 0, and that the process S 1 is the one proposed by Cartea and González-Pedraz (2012), i.e.:

S 0 (t) ≡ 0, ( 20
)
S 1 (t) = f (T ) + X(T ) + Y (T ), (21)
where S 1 (t) = S A,B (t) is the process given in equation (3).

Spread Option Pricing

In the case of a monthly interconnection product, the option gives the right to dispatch electricity from location A to location B at each hour of the considered month. Thus the valuation of interconnection capacity between two locations is equivalent to evaluating a strip of European-style options, where the underlying is the spread price between the two countries. Let r be the continuously-compounded risk-free interest rate. The price at time t of a call option on the spread S A,B at a single hour T is given by:

C A,B (t, T, K A,B ) = e -r(T -t) E t max(S A,B (T ) -K A,B , 0) .
It follows that the value of a monthly call option on the spread S A,B during the month M is then given by:

C A,B M (t, K A,B ) = TN T =T1 C A,B (t, T, K A,B ),
where T 1 ,..., T N are the hours of delivery in this month.

Using the model's characteristic function ϕ, the price of the call option can be evaluated as:

C A,B (t, T, K A,B ) = - 1 2π ∞-iα -∞-iα e -izK 1 z 2 ϕ(z)dz.
These results are well-known and stated for completeness in Appendix A.

We also have a parity between call options on spreads in opposite directions with opposite strikes. Let C A,B (t, T, K) and C B,A (t, T, -K) be two call options with the same maturity T on the opposite spreads S A,B (T ) and S B,A (T ) with opposite strikes K and -K. Their present values are then given by the expressions:

C A,B (t, T, K) = e -r(T -t) E t max(S A,B (T ) -K, 0) , C B,A (t, T, -K) = e -r(T -t) E t max(S B,A (T ) + K, 0) .
Observe that the call option C B,A (t, T, -K) with strike -K can also be seen as a put option P A,B (t, T, K) with strike K, since:

E t max(S B,A (t) + K, 0) = E t max(K -S A,B (t), 0) .
We therefore have, from standard call-put parity:

C A,B (t, T, K) -C B,A (t, T, -K) = e -r(T -t) • (E t S A,B (T ) -K).
The behaviour of electricity prices in European wholesale markets has changed significantly in recent years. Figure 4 illustrates this change, beginning around September 2021, with the example of the DE→FR spread between France and Germany. Post-covid Chinese economic revival, the European gas crisis caused by the war in Ukraine, low electricity production in France due to maintenance issues with nuclear power plants, an especially hot summer in 2021 -all of these factors have wrought profound changes in the equilibrium of commodity markets and have led to large increases in price levels and market volatility.

Our article does not attempt to provide a detailed explanation of the new fundamentals underlying this paradigm change. However, in order to estimate our model, we need to take these changes into account. We therefore split our data set into two parts: 1) from 1 January 2017 to 31 August 2021;

2) from 1 September 2021 to 31 August 2022.

From now on, we also refer to these two parts of our original data set as data set 1 and data set 2, respectively, in order to simplify our terminology.

Table 3 summarises the model parameters to be estimated. Note that we estimate them in the order given in the table, from top to bottom.

Parameter Role q 0 , q 1 infinitesimal generator Q of Markov regime-switching process M λ + , η + intensity and inverse mean of upward jumps J + λ -, η - intensity and inverse mean of downward jumps J - α, σ mean-reversion and volatility of Ornstein-Uhlenbeck process X β mean-reversion of pure jump process Y Table 3: Model Parameters.

J o u r n a l P r e -p r o o f

Journal Pre-proof

Specifying and Estimating the Seasonality Function

It is well known that electricity prices are seasonal at several frequencies, such as daily, weekly, monthly, and yearly. Accordingly, seasonality functions are an important topic in the literature; see, for example, [START_REF] Lucia | Electricity prices and power derivatives: Evidence from the Nordic power exchange[END_REF] and [START_REF] Cartea | Spot price modeling and the valuation of electricity forward contracts: The role of demand and capacity[END_REF] for discussions and various specifications.

Seasonality is less of an issue for electricity price spreads between neighbouring countries, the case considered in this article. Most of the European countries studied here are in the same time-zone and share many consumption patterns. Seasonal fluctuations therefore tend to cancel out to a large degree. And of course, as periods of market coupling increase between countries, the importance of seasonality diminishes even further, since it plays no role during coupled periods.

Moreover, the market only provides information on forward prices, which are given at a monthly granularity at the moment of the auctions. Therefore, when pricing monthly spread options, we only have a constant view of the spread for the given period. To be consistent with the pricing in the calibration procedure, we therefore shift the process by a constant whose value is simply given by the mean of the historical spread process in the data set. Since we account for regime switching of market coupling, we only consider this component during decoupled periods.

As described above, when we price spread options, we will use the market forward price that is appropriate for the option period. For example, when pricing monthly spread options, we will take the difference of the two monthly forward prices of the countries in question as input. In our model, we have:

f (T ) = π 0 f 0 (T ) + π 1 f 1 (T ).
We assume f 0 (T ) = 0 for all T . Therefore, in order for the model to reproduce the market forward f (T ) for time T , we set:

f 1 (T ) = 1 π 1 f (T ).

J o u r n a l P r e -p r o o f

Journal Pre-proof

Estimating the Transition Probabilities of the Markov Process

We continue the description of our multi-step estimation procedure by estimating the transition probabilities of the Markov process M introduced in Section 3.2.

First we need to identify the switches between the two regimes represented by M (t) = 0 and M (t) = 1. This step allows us to estimate the infinitesimal generator Q of the Markov process M . Note that this first step is directly based on the spread data and is therefore simple and stable.

Let N i→j denote the number of times the process switches from state i to state j ̸ = i, and let N i denote the number of steps spent in state i. We can then estimate the 1-hour transition probabilities as:

p 01 = p 01 (1 hour) = N 0→1 N 0 , ( 22
)
p 10 = p 10 (1 hour) = N 1→0 N 1 , (23)
which, using equations ( 17) and ( 18), gives us the estimates:

q 0 = -ln(1 -p 01 -p 10 ) p 01 p 01 + p 10 , ( 24
)
q 1 = -ln(1 -p 01 -p 10 ) p 10 p 01 + p 10 , ( 25
)
for the entries q 0 = q 01 = -q 00 and q 1 = q 10 = -q 11 of the matrix Q. Note that the coefficient 1 t = 1 in this case, since our unit of time is one hour. This estimation of the parameters of the infinitesimal generator Q completely describes the behaviour of the Markov process M over time. Table 4 shows the parameters describing the Markov process M . In the first data set, the probability of switching from the coupled regime 0 to the uncoupled regime 1 is estimated to be p 01 = 17.60% per hour. This means that the expected length of staying in the coupled regime is p -1 01 = 5.68 hours. For the second data set, this length decreases slightly to 0.1817 -1 = 5.50 hours. In the other direction, in the first data set, the probability of switching from the uncoupled regime 1 to the coupled regime 0 is 12.78% per hour, which translates into an expected length of staying in the uncoupled regime of 0.1278 -1 = 7.82 hours. For the second data set, this length increases slightly to 0.1109 -1 = 9.02 hours. Note that the transition probabilities remain relatively stable over the two data sets. We will see shortly that this is not at all the case for most of the other model parameters.

Interconnection

Figure 5 highlights the times at which regime switching occurs for the price spread between France and Germany in the first days of January 2020. The plotted spreads are in line with the random behaviour of the regime switching described in Table 4. In the first week, four longer periods of coupled markets can be observed, as well as one long period of decoupled markets. However, towards the end of the period shown, there is more frequent regime-switching, and only one slightly longer period of uncoupled markets can be seen. We apply a recursive semi-parametric filter to identify the calendar position of the jumps in the price spread. The procedure identifies the hypothetical arrival of a jump when the spread difference deviates, in absolute value, by more than three standard deviations from its mean. It is therefore dependent on the length of the estimation time-windows.

Estimating jump intensities and jump size parameters

Once we have identified all up-and down-jumps in our time series, we need to identify the respective intensities λ + and λ -of the Poisson processes N + and N -appearing in the jump process Y . We do this by determining the frequency with which the jumps occur. We then fit two exponential distributions for the size of up-and down-jumps in the jump process Y .

The price spread is highly unbalanced in the second part of our data set, with French electricity prices almost always lying above German ones. Consequently, our estimation procedure does not detect any downward jumps; downward moves of the spread are entirely due to mean-reversion. We conclude that the jump process Y only allows for upward jumps.

Our results are shown in Table 5. Regarding the first part of our data set, we see that the frequency λ + of positive jumps is higher than the frequency λ -of negative jumps. Also, the average size of positive jumps is 0.0433 -1 = 23.09 and slightly greater in magnitude than the average size of negative jumps, 0.0461 -1 = 21.69. This observation is exacerbated in the second part of our data set, where all jumps are positive, and their size is now 0.0030 -1 = 333.3 on average.

Figure 8 shows the probability density functions of the exponential distributions for the values of η + and η -given in Table 5. It can be clearly seen that the density for the second data set has a much fatter tail to the right.

J o u r n a l P r e -p r o o f

Journal Pre-proof

Estimating the Mean-Reversion Speeds and Volatility

In order to conclude our estimation of the model parameters, we need to estimate the parameters α and σ, which describe the Ornstein-Uhlenbeck process X given by equation ( 4), as well as the parameter β, which describes the mean-reversion speed of the pure jump process Y given by equation ( 5).

We proceed by minimising the mean of the squared differences between observed and modelled spreads. There is, however, a difficulty, since it is not possible to observe the realised time series followed by X and Y individually, but only their sum X + Y . We therefore first fix the mean reversion parameter β of the process Y , and then apply the least-squares method to find the parameters α and σ of the process X implied for this value of β. We then repeat this step until we find the optimal value for the parameter β, and therefore the overall optimal triplet α, σ and β.

Table 6 shows the results for the spread DE→FR for both calibration windows. The mean reversion speeds α and β are relatively stable over the two data sets. However, the estimated volatility σ of the second time window is almost 10 times higher than that of the first one. We will see in Section 5.1 below that similar, if not quite so extreme, increases in volatility also hold for other pairs of neighbouring countries. The first part of our results illustrates how dramatically electricity market conditions in Europe changed in 2021, and demonstrates that our model is capable of adapting to this new environment of high prices and volatility. The second part compares auction prices of monthly spread options from the past years to traded futures quotes and model prices for these options, and also to realised spreads during these periods.

Interconnection

A New European Paradigm

We begin with the regime-switching model estimated over data set 1, i.e., the period from 1 January 2017 to 31 August 2021. In Figure 9, we compare the realised price spread to a path that was simulated with the model. We see that our regime-switching model successfully captures the important stylised features of intermittently coupled markets:

1) two regime states, representing coupled and decoupled markets;

2) up-and downward jumps in the spread process;

3) a strong mean-reversion characteristic. We also estimate the regime-switching model over data set 2, i.e., the period from 1 September 2021 to 31 August 2022. Our estimation procedure reveals a significant increase in the spread volatility σ. Table 7 presents estimated volatility parameters for all interconnections with Germany for both data sets. We see that the volatility of the spread has been multiplied by factors ranging from close to 3 (DE → BE, DE → NL, DE → PL) to almost 10 (DE → FR). To graphically show this change in the market, Figure 10 presents three sample paths simulated for the first data set estimation on the left, and for the second data set on the right. We intentionally keep the same scale on both sides in order to highlight the radically different spread levels prevalent in the two data sets. The new paradigm in the electricity market can be clearly observed. Before the summer of 2021, as can be seen in the historical spreads on the top left, price spreads mostly stayed in the range of ±50 EUR, with a few drops to around -100 EUR. After this time, as can be seen on the top right, price spreads frequently reach levels of 300 EUR, and the tenfold increase in volatility is easily visible. Furthermore, the simulations show that our model is able to capture the unbalanced spread that we observe in the historical data. Note that although the estimated downward jump intensity λ -is equal to zero for the second data set, we still observe negative spreads in the simulations on the right-hand side of Figure 10. This is because the continuous part of the spread process, X, is not defined individually for negative and positive spreads, and can therefore in itself lead to negative spreads.

A Comparison of Traded and Real Options

Table 8 compares market and model prices of monthly spread options for the country couple (DE, FR) in both directions. The market prices were obtained from the Joint Allocation Office (JAO) "Auctions" website. We also report monthly futures prices from the European Energy Exchange (EEX) for the closing date of the JAO auction for reference, together with their difference or spread. It is important to note that these spread option prices are auction prices, obtained from a 3-day bidding period and published on the closing date. The link between option prices and futures prices is therefore less direct here than in situations where option prices are real-time trading prices, and futures and options can be traded in synchrony.

Tables 9 and10 show the profit-and-loss results of a trading strategy that simply makes bids for monthly spread options in the JAO auction at the model price. We assume here for simplicity that our bid does not alter the final auction price, which is reasonable if the volume is low compared to the total volume of interconnection exchange. Table 9 shows the results for the model estimated on the first data set, and Table 10 for the model estimated on the second data set. Naturally, in order to obtain proper out-of-sample results, we only use a model for bidding at auctions occurring after the end of the model's estimation period. If our bid succeeds in the auction, i.e., our model price is higher than the auction price, then we make a profit if the spread that is realised during the course of the calendar month is higher than the auction price at which we bought the option; otherwise, we make a loss. If our bid is too low and therefore not executed in the auction, we still compare the auction price to the realised spread and report whether our insufficient bid turned out to be right or wrong.

The analysis of our trading results is revealing and sheds further light on our findings about the "new European paradigm" in Section 5.1. Table 9 shows that the trading strategy executed with the model estimated for data set 1 is highly successful from September 2021 to August 2022: 20 out of 24 bids are made in the correct direction, and the overall profit is decidedly positive. However, with the market turmoil beginning in the summer of 2022, the situation reverses completely, and the trading strategy leads to injudicious bids and significant losses from September 2022 onwards. Interestingly, the trading outcome is essentially the same when we use prices from the model estimated on the more recent data set 2: bids are now also made for the direction FR → DE and lead to a slight loss; however, the bids in the more expensive direction DE → FR are unchanged and lead to the same losses as with the first model.

The finding that our model, used with two starkly differing sets of parameters, leads to essentially identical trading results is in part explained by the futures spreads of Table 8. These futures spreads turn out to be much higher than the realised spreads from September to December 2022. And since the futures spreads are a first-order input to the pricing model, it is a direct consequence that the model prices are also much higher than the realised spreads.

Finally, we can see from Tables 9 and 10 that from February 2023 onward, JAO auction prices have come down and realigned with the realised spreads. From February to August 2023, Model 1 performs marginally worse than Model 2 in the direction DE → FR (-6.95 vs -6.75 EUR), and slightly better (2.47 vs 1.98 EUR) in the direction FR → DE. Given that the realised spreads have decreased substantially in this period, we believe that it is necessary to regularly re-estimate model parameters on a rolling window basis. This allows us to work with a model that is a good fit with current market conditions and that can achieve positive trading results, as for example with Model 1 in the period before September 2022. Auctions of transmission rights between neighbouring countries are becoming increasingly active. In a parallel development, the introduction of market coupling frequently leads to smaller price differences between such countries. Indeed, if two countries are completely coupled, the price of a given hour of electricity will be identical in each country, resulting in a price spread of zero.

Clearly, it is important to take this market coupling into account when evaluating transmission rights, as neglecting it would lead to a significant overvaluation of this right. In order to address this issue, we introduce a general regime-switching mechanism that can be applied to many models in the literature. We describe the model estimation procedure in detail, and compare model and market prices of European spread options. We observe a dramatic paradigm shift in our data set at the end of the summer of 2021, and show that this shift has a strong effect on the model parameters. We also see that the reliable pricing and trading of spread options has become more problematic.

In future work, it would be interesting to investigate the causes of regime-switching, and introduce the possibility of correlating the states of the regime-switching Markov process to electricity price levels or to other variables. These variables could include market volatility, prices in related markets such as natural gas markets and weather conditions such as sun and wind, which affect the electricity produced by renewable sources.

A The Call Option Formula via the Fourier Transform

The characteristic function ϕ = ϕ T of the spot price S T at time T is given for u ∈ R by:

ϕ T (u) = E[e iuST ] = ∞ -∞ e iuS p(S)dS, ( 26
)
where p is the probability density function of S T .

For example, for the Bachelier model:

dS t = µdt + σdB t , S 0 ∈ R,
with drift µ ∈ R and volatility σ > 0, and with solution:

S T = S 0 + µT + σB T , B T ∼ N (0, T ), the characteristic function ϕ B is given by: ϕ B (u) = e iu(S0+µT )-1 2 u 2 σ 2 T .

In general, for a European call option with strike K and maturity T , we have:

C(K) = e -rT ∞ K (S -K)p(S)dS.
For α > 0, define the modified or damped call price by:

c(K) = e αK C(K). (27)

Figure 1 :

 1 Figure 1: Price differences between France and Germany since 1 January 2017.

Figure 2 :

 2 Figure 2: Hourly price differences between France and Germany in May 2018.

Figure 3 :

 3 Figure 3: Daily price differences between France and Germany in the first half of 2018. The daily spread was only equal to zero on 16 April 2018.

ϕ

  T (u) = E[e iuS(T ) ] = P(M T = 0)E[e iuS(T ) |M T = 0] + P(M T = 1)E[e iuS(T ) |M T = 1] = P(M T = 0)E[e iuS0(T ) ] + P(M T = 1)E[e iuS1(T ) ]

Figure 4 :

 4 Figure 4: Spread DE→FR before and after September 2021.

Figure 6 :

 6 Figure 6: Process netted of points from regime switching.

Figure 7

 7 Figure 7 highlights the jumps we identified for the DE→FR spread in the first half of 2020, during the estimation procedure for the first part of our data set (from 01/01/2017 to 31/08/2021).

Figure 8 :

 8 Figure 8: Exponential densities for both data sets.

Figure 9 :

 9 Figure 9: Comparison of realised and simulated price spreads for data set 1.

Table 1 :

 1 German (DE) price spreads with neighbouring countries for 24 May 2022 HOUR DE → AT DE → BE DE → FR DE → NL DE → PL DE → CZ

	20.97	69.21	80.93	137.93	16.34	16.34
	21.25	48.02	82.91	93.75	21.25	21.25
	23.61	42.88	72.68	75.93	15.55	15.55
	17.46	33.28	48.64	41.76	0	0
	8.35	20.2	29.64	25.43	0	0
	0	0	0	0	-37.58	0
	0	0	0	0	-60.36	0
	0	0	0	0	-95.82	0
	0.66	-2.04	2.86	-1.88	-89.47	0.66
	0	0	0	0	-70.51	0
	4.6	-3.11	0.13	-3.7	-60.27	4.6
	16.64	0	6.11	0	-32.75	16.64
	19.54	0	50.83	0	-1.44	19.54
	32.45	0	50.69	0	11.92	32.45
	35.5	0	50.9	0	32.7	35.5
	15.86	0	10.95	0	15.86	15.86
	7.6	-4.28	0	-4.28	7.6	7.6
	13.57	-7.65	-1.71	-6.19	23.52	23.52
	0	0	0	0	10.83	0
	0	0	0	0	30.45	0
	0	0	0	0	0	0
	0	0	0	0	0	0
	0	0	0	0	0	0
	0	0	0	0	-15.12	0

Table 2 :

 2

	Peak	146.42	0	82.91	0	0	156.37
	Off-Peak	91.64	213.41	314.80	374.80	0	53.14
	Base	238.06	196.33	485.56	358.75	0	209.51
	Hourly	238.06	213.41	487.27	374.80	186.02	209.51

Interconnection option payoffs for 24 May 2022 Method DE → AT DE → BE DE → FR DE → NL DE → PL DE → CZ

Table 4 :

 4 Transition probabilities of the Markov process M .

		DE→FR
	Start Date	01/01/2017 01/09/2021
	End Date	31/08/2021 31/08/2022
	p 01 (1 hour)	0.1760	0.1817
	p 10 (1 hour)	0.1278	0.1109
	q 0	0.2098	0.2150
	q 1	0.1523	0.1312

Table 6 :

 6 Mean reversion and volatility parameters of the processes X and Y .

		DE→FR
	Start Date	01/01/2017 01/09/2021
	End Date	31/08/2021 31/08/2022
	α	0.060	0.083
	σ	4.506	43.632
	β	0.100	0.100

Table 7 :

 7 German (DE) spread volatility σ with neighbouring countries for both data sets. data set DE → AT DE → BE DE → FR DE → NL DE → PL DE → CZ

	1	3.022	6.543	4.506	3.676	11.584	3.010
	2	14.060	17.306	43.632	11.021	33.980	14.283

Table 8 :

 8 Market and Regime-Switching Model Prices of Monthly DE-FR Spread Options in EUR/MWh

	JAO	Market			JAO
	Auction	Futures			Options		Model 1 Model 2 Model 1 Model 2
	Closing Date	DE	FR Spread DE→FR FR→DE DE→FR DE→FR FR→DE FR→DE
	22/09/21	135.07 154.26	19.19	18.27	1.08	20.44	-	0.21	-
	25/10/21	191.42 228.80	37.38	37.06	1.29	38.59	-	0.03	-
	24/11/21	186.62 268.23	81.61	66.13	1.03	82.82	-	0.00	-
	23/12/21	346.28 431.23	84.95	40.00	1.31	86.16	-	0.00	-
	24/01/22	210.97 251.38	40.41	39.75	0.85	41.62	-	0.03	-
	23/02/22	198.83 225.93	27.10	29.11	1.01	28.32	-	0.09	-
	23/03/22	242.39 260.00	17.61	52.08	0.85	18.88	-	0.26	-
	22/04/22	205.14 226.23	21.09	21.00	1.59	22.33	-	0.17	-
	24/05/22	184.65 208.55	23.90	23.03	0.87	25.13	-	0.12	-
	22/06/22	274.51 316.97	42.46	44.10	0.80	43.67	-	0.02	-
	25/07/22	351.42 411.59	60.17	63.08	0.62	61.38	-	0.00	-
	23/08/22	550.56 650.00	99.44	95.78	2.41	100.65	-	0.00	-
	22/09/22	365.61 475.87 110.26	100.00	1.13	111.47	116.61	0.00	3.73
	24/10/22	243.44 435.33 191.89	180.00	1.55	193.10	196.66	0.00	1.30
	23/11/22	276.98 442.50 165.52	178.60	1.33	166.73	170.47	0.00	1.68
	15/12/22	327.07 416.87	89.80	91.54	1.38	91.01	97.58	0.00	5.51
	24/01/23	131.50 152.12	20.62	19.56	1.16	21.86	42.94	0.18	21.72
	23/02/23	129.61 146.76	17.15	15.11	0.93	18.42	40.78	0.28	23.13
	23/03/23	106.52 115.64	9.12	8.98	1.55	10.77	36.04	0.97	26.65
	25/04/23	97.62	99.17	1.55	4.65	2.12	4.93	31.91	3.24	30.32
	23/05/23	83.65	80.17	-3.48	2.61	5.49	2.43	29.35	6.21	32.93
	23/06/23	95.07	96.96	1.89	3.59	3.67	5.14	32.09	3.09	30.14
	25/07/23	89.44	81.77	-7.67	2.77	9.15	1.23	27.33	9.50	35.22

Table 9 :

 9 Performance of trading strategy calculated with model for data set 1

	JAO	Realised		JAO				Trading
	Auction	Spread		Price		Model 1		Strategy
	Month DE→FR FR→DE DE→FR FR→DE DE→FR FR→DE DE→FR FR→DE
	Oct 21	33.55	0.67	18.27	1.08	20.44	0.21	15.28 not exec
	Nov 21	41.58	0.67	37.06	1.29	38.59	0.03	4.52 not exec
	Dec 21	53.99	0.37	66.13	1.03	82.82	0.00	-12.14 not exec
	Jan 22	43.87	0.17	40.00	1.31	86.16	0.00	3.87 not exec
	Feb 22	56.75	0.01	39.75	0.85	41.62	0.03	17.00 not exec
	Mar 22	43.69	0.49	29.11	1.01	28.32	0.09 not exec not exec
	Apr 22	67.50	0.13	52.08	0.85	18.88	0.26 not exec not exec
	May 22	20.12	0.17	21.00	1.59	22.33	0.17	-0.88 not exec
	Jun 22	31.06	0.69	23.03	0.87	25.13	0.12	8.03 not exec
	Jul 22	85.90	0.04	44.10	0.80	43.67	0.02 not exec not exec
	Aug 22	28.46	1.16	63.08	0.62	61.38	0.00 not exec not exec
	Sep 22	49.14	0.55	95.78	2.41	100.65	0.00	-46.64 not exec
	Oct 22	28.22	2.02	100.00	1.13	111.47	0.00	-71.78 not exec
	Nov 22	19.15	0.90	180.00	1.55	193.10	0.00	-160.85 not exec
	Dec 22	20.07	0.80	178.60	1.33	166.73	0.00 not exec not exec
	Jan 23	14.64	0.37	91.54	1.38	91.01	0.00 not exec not exec
	Feb 23	20.61	0.16	19.56	1.16	21.86	0.18	1.05 not exec
	Mar 23	10.38	0.94	15.11	0.93	18.42	0.28	-4.73 not exec
	Apr 23	6.20	0.58	8.98	1.55	10.77	0.97	-2.78 not exec
	May 23	2.57	6.74	4.65	2.12	4.93	3.24	-2.08	4.62
	Jun 23	1.18	4.64	2.61	5.49	2.43	6.21 not exec	-0.85
	Jul 23	5.18	5.14	3.59	3.67	5.14	3.09	1.59 not exec
	Aug 23	4.40	7.85	2.77	9.15	1.23	9.50 not exec	-1.30

Table 10 :

 10 Performance of trading strategy calculated with model for data set 2

	JAO	Realised		JAO				Trading
	Auction	Spread		Price		Model 2		Strategy
	Month DE→FR FR→DE DE→FR FR→DE DE→FR FR→DE DE→FR FR→DE
	Oct 22	28.22	2.02	100.00	1.13	116.61	3.73	-71.78	0.89
	Nov 22	19.15	0.90	180.00	1.55	196.66	1.30	-160.85 not exec
	Dec 22	20.07	0.80	178.60	1.33	170.47	1.68 not exec	-0.53
	Jan 23	14.64	0.37	91.54	1.38	97.58	5.51	-76.90	-1.01
	Feb 23	20.61	0.16	19.56	1.16	42.94	21.72	1.05	-1.00
	Mar 23	10.38	0.94	15.11	0.93	40.78	23.13	-4.73	0.01
	Apr 23	6.20	0.58	8.98	1.55	36.04	26.65	-2.78	-0.97
	May 23	2.57	6.74	4.65	2.12	31.91	30.32	-2.08	4.62
	Jun 23	1.18	4.64	2.61	5.49	29.35	32.93	-1.43	-0.85
	Jul 23	5.18	5.14	3.59	3.67	32.09	30.14	1.59	1.47
	Aug 23	4.40	7.85	2.77	9.15	27.33	35.22	1.63	-1.30

This situation is similar for other European interconnections. Table1shows German (DE) hourly

J o u r n a l P r e -p r o o f

Journal Pre-proof

Detecting Jumps and Estimating their Intensity and Distribution

The goal of this section is to estimate the parameters λ + , η + and λ -, η -, which describe the intensity and inverse mean of the upward and downward jumps J + and J -given in equations ( 6) and ( 7), respectively. The remaining mean-reversion parameter β, needed to completely describe the pure jump process Y given in equation (5), will be estimated in the following Section 4.5.

In the previous Section 4.3, we identified the state of the Markov process M at all times (hours) t in our sample. We remove all t for which M (t) = 0, i.e., S(t) = S 0 (t) = 0, from our spread series S(t) in order to focus on the estimation of the parameters of the process S 1 from now on.

In this new, reduced series, however, we want to avoid counting spurious jumps that result from gluing the series together after having taken out an interval during which S(t) = 0. For example, let us assume that the last time the series is in state 1 before a regime switch 1 → 0 is t i , and then again the first time the series is in state 1 after a regime switch 0 → 1 is t j ; then, after removing the data points for times t i+1 , ..., t j-1 from our series, we tag the price move from S(t i ) to S(t j ) as an invalid candidate for a jump.

Figure 6 illustrates the process with these points removed. In the upper graph, we have marked in red the "false" jumps, caused by regime-switching, that we wish to exclude from the estimation of the jump process Y . In the lower graph, we show the sequence of price spreads for which we are actually going to estimate the jump parameters λ + , η + and λ -, η -.

Detecting potential jumps

To detect big price movements and potential jump-times, and then determine whether such a big movement is really due to a jump or whether it is only an effect of mean-reversion, we follow the methodology presented in Cartea and González-Pedraz (2012).

J o u r n a l P r e -p r o o f

Journal Pre-proof

J o u r n a l P r e -p r o o f

Journal Pre-proof

Then the Fourier transform ĉ of c can be calculated as:

The inner integral is calculated using integration by parts as:

Using this result, we finally obtain:

The call option price is then obtained via Fourier inversion as:

Since the imaginary part of the integrand is odd in u, and the real part even, this integral can also be calculated as:

Note that the function ĉ has a pole of order 2 at u = iα.

An equivalent formulation is often found in the literature, e.g., [START_REF] Lee | Option pricing by transform methods: Extensions, unification, and error control[END_REF] and Cartea and González-Pedraz ( 2012), which we will present next. By introducing the variable z = uiα, and using iz = iu + α, we can define the function Ĉ as:

J o u r n a l P r e -p r o o f

Journal Pre-proof

In terms of the transform Ĉ, the call option price is expressed as:

We denote the integrand above as a function g, g(z) := e -izK 1 z 2 ϕ(z), and proceed to calculate the residue of g at z = 0, or u = iα in terms of the variable u. Using (26) together with the exponential series, we can expand g as:


Youssouf Mini

Ky Ngac

Eric Nguyen

Semail

Sensorless Control for Non-Sinusoidal Five-Phase Interior PMSM Based on Sliding Mode Observer

Keywords: back-EMF observer, electrical integrated drive, five-phase interior permanent magnet synchronous machine, sensorless control, sliding mode observer

This paper proposes a sensorless control strategy based on Sliding Mode Observer (SMO) for a Five-phase Interior Permanent Magnet Synchronous Machine (FIPMSM), with a consideration of the third harmonic component. Compared to conventional three-phase machines, the third harmonic of back electromotive force (back-EMF) contains more information. Thus, in this paper, the first and third harmonic components of the five-phase machine are considered to estimate the rotor position which is necessary for the vector control. Simulation results are shown to verify the feasibility and the robustness of the proposed sensorless control strategy.

INTRODUCTION

Multiphase machines present some advantages compared to conventional three-phase machines, such as compactness, reliability (operating under the loss of one or more phases), and a reduction in torque ripple at low frequencies even with non-sinusoidal back-EMFs [START_REF] Levi | Multiphase electric machines for variable-speed applications[END_REF][START_REF] Khan | On performance figures of multiphase machines[END_REF]. Recently, multiphase machines have been used in electric drives with the power inverters integrated in machines [START_REF]Intelligent Integrated Energy Converter project (CE2I)[END_REF][START_REF] Burkhardt | Design and analysis of a highly integrated 9-phase drivetrain for EV applications[END_REF]. The main advantage of this integration is to reduce the global volume and weight of the integrated drives without electromagnetic compatibility phenomena [START_REF] Jahns | The past, present, and future of power electronics integration technology in motor drives[END_REF][START_REF] Abebe | Integrated motor drives: State of the art and future trends[END_REF]. In fact, this can be an effective solution for applications which require high power density and compactness, such as automotive, marine and aerospace applications [START_REF] Hansen | History and state of the art in commercial electric ship propulsion, integrated power systems, and future trends[END_REF].

In this context, the replacement of the position encoder mounted at the end of the rotor shaft by a soft position sensor using only already integrated electrical or magnetic sensors becomes interesting. Besides, to improve the precision, which is classically the weakness of the soft sensor, an algorithm taking advantage of the specificities of non-sinusoidal multiphase machines is necessary. In addition, the multiphase machines are appreciated for their tolerance, a soft position sensor added to the position encoder can be also required to bestow redundancy for the angular position used in the vector control.

In the literature, several studies have proposed sensorless control methods for interior permanent magnet synchronous machines (IPMSMs) [START_REF] Vas | Sensorless Vector and Direct Torque Control[END_REF][START_REF] Benjak | Review of position estimation methods for IPMSM drives without a position sensor part I: Nonadaptive methods[END_REF]. Many sensorless control studies are based on observer for three-phase IPMSMs, however only few papers have considered the sensorless control for multiphase machines [START_REF] Zhang | Fault-tolerant sensorless control of a five-phase FTFSCW-IPM motor based on a wide-speed strongrobustness sliding mode observer[END_REF][START_REF] Ramezani | The modeling and position-sensorless estimation technique for a nine-phase interior permanent-magnet machine using highfrequency injections[END_REF][START_REF] Parsa | Sensorless direct torque control of five-phase interior permanent-magnet motor drives[END_REF]. Several methods based on the observer, as Model Reference Adaptive System (MRAS) [START_REF] Sun | A 2-D fuzzy logic based MRAS scheme for sensorless control of interior permanent magnet synchronous motor drives with cyclic fluctuating loads[END_REF], Extended Kalman Filter (EKF) [START_REF] Michalski | Sensorless control of five phase PMSM based on extended Kalman filter[END_REF], Luenberger Observer (LO) [START_REF] Zhang | A comparative study of pmsm sensorless control algorithms: Model based vs luenberger observer[END_REF] and Sliding Mode Observer (SMO) [START_REF] Hosseyni | New sensorless sliding mode control of a five-phase permanent magnet synchronous motor drive based on sliding mode observer[END_REF] can be used to perform the sensorless control of electrical machines. Among these methods, SMO will be chosen to achieve the sensorless control of FIPMSMs due to its simple implementation compared to EKF, which results in the calculation burden especially in the case of multiphase machines. Furthermore, in terms of robustness, SMO presents a robust structure against variations of machine parameters (that will be highlighted in this paper) and noise compared to MRAS and LO [START_REF] Hosseyni | Sensorless sliding mode observer for a five-phase permanent magnet synchronous motor drive[END_REF][START_REF] Qiao | New sliding-mode observer for position sensorless control of permanent-magnet synchronous motor[END_REF]. Recently, several researches have proved that the chattering phenomenon inherent in the SMO (main disadvantage of SMO who cannot be completely eliminated) can be reduced, by replacing the saturation function by a sigmoid function [START_REF] Qiao | New sliding-mode observer for position sensorless control of permanent-magnet synchronous motor[END_REF].

In non-sinusoidal multiphase machines, the torque is produced by several harmonics of their back-EMFs and currents. The current regulation requires the rotor position to perform the vector control. Therefore, the estimation of rotor position through each harmonic (that produces a torque) can increase the degree of freedom for the current control loop, and it allows to separately control different harmonics. This can improve the reliability of the control system of FIPMSMs.

Therefore, in this paper, an observer based on sliding mode will be implemented to achieve the sensorless control of considered FIPMSM. The main contribution of this paper is to use not only one harmonic of the FIPMSM as in [START_REF] Zhang | Fault-tolerant sensorless control of a five-phase FTFSCW-IPM motor based on a wide-speed strongrobustness sliding mode observer[END_REF][START_REF] Hosseyni | Sensorless sliding mode observer for a five-phase permanent magnet synchronous motor drive[END_REF][START_REF] Yang | Iterative sliding mode observer for sensorless control of five-phase permanent magnet synchronous motor[END_REF], but also two characteristic harmonics of the five-phase machine to achieve the sensorless control. Therefore, the first and third harmonic components of the back-EMF of the FIPMSM are used to estimate the rotor position, to perform an accurate vector control.

MATHEMATICAL MODEL OF FIVE-PHASE IPMSM

The FIPMSM model in natural frame, without the magnetic saturation and saliency, is given by [START_REF] Semail | Right harmonic spectrum for the back-electromotive force of an n-phase synchronous motor[END_REF]:

  di v Ri e dt = + + L (1)
with:

1 2 3 4 5 T v v v v v v =   1 2 3 4 5
T e e e e e e =  

1 2 3 4 5
T i i i i i i =         =             T ( 2
)
where, 𝛼 = 2𝜋

5

. The model of the FIPMSM in the stationary reference frame (α-β) is given by the following equation:

           + - - = + - - = + - - = + - - = 3 3 3 3 3 3 3 3 1 1 1 1 1 1 1 1                 v e Ri dt di L v e Ri dt di L v e Ri dt di L v e Ri dt di L s s p p (3)
where, 𝑖 ⃗ 𝛼𝛽1 = [𝑖 𝛼1 𝑖 𝛽1 ] 𝑇 and 𝑖 ⃗ 𝛼𝛽3 = [𝑖 𝛼3 𝑖 𝛽3 ] 𝑇 represent respectively currents of main and secondary fictitious machines. Lp and Ls represent respectively inductances of main and secondary fictitious machines. As the fictitious machines are mechanically coupled, the electromagnetic torque of the FIPMSM can be obtained by the sum of torques provided by all fictitious machines. Therefore, the total torque is calculated as follows: . where, Γ1 is the torque of the main fictitious machine; Γ3 is the torque of the second fictitious machine; and Ω is the mechanical speed.

The back-EMF in stationary reference frame (α-β) of each fictitious machine is expressed as (only the 1 st and 3 rd harmonics are considered) [START_REF] Semail | Vectorial formalism for analysis and design of polyphase synchronous machines[END_REF]:

            =-   =   =-   =-  (5)
where, ψ1 and ψ3 are respectively the first and third harmonic components of permanent magnet flux linkage, ωr is the electrical angular velocity. The angles θp and θs are defined as: θp=pΩt+θ0p and θs=3pΩt+θ0s. Where θ0p and θ0s represent respectively the initial angles of the first and third harmonics of the back-EMF.

To perform accurate vector control, the rotor position and speed information are required to compute the Park transformation in the rotor reference frame. In this context, it can be seen in ( 5) that the back-EMF signal contains this two information. From [START_REF] Jahns | The past, present, and future of power electronics integration technology in motor drives[END_REF], it can be noticed that in the case of the non-sinusoidal FIPMSM (where the back-EMF contains the 1 st and 3 rd harmonic), the rotor position and speed can be estimated through the two fictitious machines (defined by the Concordia matrix (2)). Thus, an observer based on Sliding Mode will be designed, in order to estimate with high accuracy, the back-EMF signals necessary to extract the rotor position and speed information. It should be noted that the angles required to control the fictitious machines can be estimated by: using only the estimated back-EMF signals of 1 st harmonic, or all estimated back-EMF signals of 1 st and 3 rd harmonics. The two approaches will be discussed in section 3.3.

SLIDING MODE OBSERVER DESIGN

Current observer

The Sliding Mode Observer (SMO) can be designed in the stationary reference frame (α-β) [START_REF] Kim | A high-speed slidingmode observer for the sensorless speed control of a PMSM[END_REF]. It is based on the measured stator currents, and the calculated (estimated) ones by the mathematical model of the machine. The SMO is constructed by comparing the measured stator current at the estimated stator current in (α-β) frame. Indeed, the purpose is to minimize the error between the measured and the estimated stator current by using a switching function (as saturation function, sign function or sigmoid function) [START_REF] Chi | Implementation of a sliding-mode-based position sensorless drive for highspeed micro permanent-magnet synchronous motors[END_REF].

We define the error vector 𝑆 ⃗ , which belongs to the sliding surface when 𝑖⃗ 𝑠 ≃ 𝑖 ⃗ 𝑠 , as:

ˆSS S i i =- (6)
where, 𝑖 𝑠 ⃗⃗⃗ = [𝑖 𝛼1 𝑖 𝛽1 𝑖 𝛼3 𝑖 𝛽3 ] 𝑇 is the estimated stator current vector and 𝑖 𝑠 ⃗⃗⃗ = [𝑖 𝛼1 𝑖 𝛽1 𝑖 𝛼3 𝑖 𝛽3 ] 𝑇 is the actual measure of the stator current vector. the vector 𝑆 ⃗ is defined for the FIPMSM as:

1 11 1 11 3 33 3 33 Ŝ ii S ii S S ii S ii           -   -   ==   -   -     (7)
By using the mathematical model of the FIPMSM (3) and the sliding mode theory, the current observer based on SMO can be designed as follows:

1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 ˆp p s s di L Ri v z dt di L Ri v z dt di L Ri v z dt di L Ri v z dt                   = - + -        = - + -          = - + -        = - + -      (8)
where we assume that Zα1, Zβ1, Zα3 and Zβ3 represent the outputs of switching functions that contain the back-EMF signal and a high frequency component.

( ) (

)

1 1 1 1 1 1 1 1 ẑ k F i i z k F i i        =-   =-   ( ) (
)

3 2 3 3 3 2 3 3 ẑ k F i i z k F i i        =-   =-  
k1 and k2 represent the constant observer gains. The saturation function or sign function used in the conventional SMO are replaced by a continuous function, i.e., the sigmoid function, which is defined as:

𝐹(𝑥) = [ 2 (1+𝑒 -𝑎𝑥 ) ] -1.
Where x is a variable, and a is the positive adjustable parameter for the slope of the sigmoid function.

The aforementioned observer of current based on SMO is stable if it converges toward the sliding surface, where the error is equal to zero. Indeed, to verify the stability of the current observer, the Lyapunov function is utilized.

The Lyapunov function is selected as:

1 1 3 3 1 1 3 3 1 2 T V S S S S S S S S             =     (9)
Two conditions are required to guarantee the stability of the aforementioned SMO:

✓ The Lyapunov function is positive definite. ✓ The derivative of Lyapunov function should be negative.

0 V  (10)
The error equation is obtained by subtracting (3) from ( 8) as follows:

( ) ( ) ( ) ( ) 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 2 3 3 3 3 3 2 3 3 ˆˆˆp p s s dS L RS e k F i i dt dS L RS e k F i i dt dS L RS e k F i i dt dS L RS e k F i i dt                       = - + - -       = - + - -        = - + - -       = - + - -     (11)
The Lyapunov function V is positive definite, because it's the sum of the square of the stator current in (α-β) frame [START_REF] Yang | Iterative sliding mode observer for sensorless control of five-phase permanent magnet synchronous motor[END_REF]. Therefore, to guarantee the stability condition of SMO, it is only needed to prove that the derivative of Lyapunov function is negative. It can be expressed as:

1 1 3 3 1 1 3 3 1 1 1 1 3 3 3 3 0 T V S S S S S S S S S S S S S S S S                   =   = + + + 
The condition in ( 10) is satisfied if k1 and k2 are large enough, and bounded as in [START_REF] Hosseyni | Sensorless sliding mode observer for a five-phase permanent magnet synchronous motor drive[END_REF][START_REF] Qiao | New sliding-mode observer for position sensorless control of permanent-magnet synchronous motor[END_REF][START_REF] Yang | Iterative sliding mode observer for sensorless control of five-phase permanent magnet synchronous motor[END_REF]:

( ) (
)

1 1 1 2 3 3 max , max , k e e k e e          (12)
So, the stability of the current observer based on the sliding mode observer will be guaranteed by choosing the appropriate gains k1 and k2.

Back-EMF observer

Based on the current observer (8), the equivalent back-EMF signal can be obtained through the output of the switching function. But the signal still contains high frequency components and cannot be used for the estimation of the rotor position and speed. Thus, to extract the back-EMF signals, an observer will be elaborated [START_REF] Hosseyni | Sensorless sliding mode observer for a five-phase permanent magnet synchronous motor drive[END_REF][START_REF] Qiao | New sliding-mode observer for position sensorless control of permanent-magnet synchronous motor[END_REF][START_REF] Yang | Iterative sliding mode observer for sensorless control of five-phase permanent magnet synchronous motor[END_REF].

It is assumed that the speed changes slowly (the derivative of the rotor electrical angular velocity can be considered equal to zero approximatively 𝜔̇𝑟 = 0). Based on (5), the back-EMF model of each fictitious machine can be expressed as:

                         = - = -    = - =     = - =    = = -   (13)
From ( 13) and the current observer in [START_REF] Vas | Sensorless Vector and Direct Torque Control[END_REF], the back-EMF observer for the sensorless control of the FIPMSM is written as follows:

( ) ( ) ( ) ( ) 1 1 1 1 1 1 1 1 1 1 3 3 2 3
                     = - - -    = - -     = - -    = - - -   (14)
where, 𝑒̂𝛼 1 , 𝑒̂𝛽 1 , 𝑒̂𝛼 3 and 𝑒̂𝛽 3 are the estimated values of the back-EMF in the stationary reference frame (α-β), l1 and l2 are the constant gains which are determined through the stability conditions according to the Lyapunov function, in the same way of the aforementioned current observer.

Based on the Ref. [START_REF] Hosseyni | Sensorless sliding mode observer for a five-phase permanent magnet synchronous motor drive[END_REF][START_REF] Qiao | New sliding-mode observer for position sensorless control of permanent-magnet synchronous motor[END_REF][START_REF] Yang | Iterative sliding mode observer for sensorless control of five-phase permanent magnet synchronous motor[END_REF], the observer gains values l1 and l2 should be greater than zero to guarantee the stability of the back-EMF observer.

Rotor position and speed estimate

The rotor position and speed, required to achieve the accurate sensorless control of the FIPMSM, are estimated through the extracted back-EMF signals (that contain the first and third harmonic). Therefore, by using the back-EMF estimated from the observer based on sliding mode [START_REF] Michalski | Sensorless control of five phase PMSM based on extended Kalman filter[END_REF], and the relationship between the back-EMF and the rotor position as shown in [START_REF] Jahns | The past, present, and future of power electronics integration technology in motor drives[END_REF], the estimation value of the rotor position is given in [START_REF] Hosseyni | New sensorless sliding mode control of a five-phase permanent magnet synchronous motor drive based on sliding mode observer[END_REF]. From ( 5), the electrical rotor speed can be estimated using the estimated back-EMF of 1 st harmonic as:

 =-    =     =-    =-   (16)
Conventional SMO, developed in [START_REF] Zhang | Fault-tolerant sensorless control of a five-phase FTFSCW-IPM motor based on a wide-speed strongrobustness sliding mode observer[END_REF][START_REF] Hosseyni | Sensorless sliding mode observer for a five-phase permanent magnet synchronous motor drive[END_REF][START_REF] Qiao | New sliding-mode observer for position sensorless control of permanent-magnet synchronous motor[END_REF][START_REF] Yang | Iterative sliding mode observer for sensorless control of five-phase permanent magnet synchronous motor[END_REF], uses only the 1 st harmonic component to compute (𝑠𝑖𝑛𝜃 ̂𝑝, 𝑐𝑜𝑠𝜃 ̂𝑝) and then (𝑠𝑖𝑛𝜃 ̂𝑠, 𝑐𝑜𝑠𝜃 ̂𝑠), in order to control the main and secondary fictitious machine of the FIPMSM. In fact, angle 𝜃 ̂𝑠 is computed by multiplying 𝜃 ̂𝑝 by 3. However, in [START_REF] Jahns | The past, present, and future of power electronics integration technology in motor drives[END_REF], angle θp contains an offset θ0p which could be different from an offset θ0s contained in θs. This means that multiplying 𝜃 ̂𝑝 by 3 does not allow to obtain the real angle 𝜃 ̂𝑠. To avoid this problem, each angle should be estimated from the corresponding back-EMF signals. For this purpose, the proposed sensorless control uses the 1 st harmonic component to estimate only (𝑠𝑖𝑛𝜃 ̂𝑝, 𝑐𝑜𝑠𝜃 ̂𝑝) , and the 3 rd harmonic to estimate (𝑠𝑖𝑛𝜃 ̂𝑠, 𝑐𝑜𝑠𝜃 ̂𝑠). This approach can guarantee the independence between the control loops of different fictitious machines.

It should be noted that as the back-EMF of 1 st harmonic is more important than the 3 rd harmonic one, the estimation process of rotor speed by 1 st harmonic is more accurate. The overall block diagram of sensorless control of the FIPMSM is shown in Figure 1. The estimated rotor position, by using all fictitious machines, can also be expressed as:

SIMULATIONS RESULTS

Based on the aforementioned sliding mode observer, the simulation is used to verify the effectiveness of the proposed SMO for the sensorless control of the FIPMSM. The proposed system from Figure 1 has been implemented in the MATLAB/Simulink programming environment. The PWM switching frequency is 10kHz. The sampling time used for the sensorless control system shown in Figure 1 is set at 1 µs. It can be noted that the low and zero speed region is not considered in this study. In this simulation, a special cycle for the reference rotor speed as shown in Figure 2 (a) is considered to verify the stability and robustness of the proposed SMO under the speed and load torque variations. The FIPMSM parameters are provided in Table 2. The SMO parameters are provided in Table 3. The control strategy with id=0 is carried out. In Figure 2 In Figures 3 and4, the current and back-EMF observers based on SMO accurately estimate the current and the back-EMF signals of each fictitious machine in wide speed range. Thus, the stability and robustness of the proposed SMO under the torque and speed variations are proved. When the amplitude of the back-EMF signal of the main and secondary fictitious machine is low, the estimation process is not precise. Therefore, the estimation of the rotor position in this range will be considered greatly impacted.

Feasibility of the proposed sensorless control

In Figure 2 (a), the reference rotor speed cycle contains several transient and steady states, which are used to highlight the effectiveness of the proposed SMO. The estimated currents and back-EMF signals show that the proposed SMO is not impacted when the reference speed changes from steady state to transient state or vice versa.

It is noted that the load torque disturbance also has no obvious effect on the estimation process. The load torque is applied in transient state as in steady state. In fact, it can be seen form the simulation waveforms of Figure 4 that the estimated back-EMF signals of the first and third harmonic components are not affected by the load torque variations.

In Figure 3 (c), the current of the secondary fictitious machine is less than the current in the main fictitious machine. Therefore, the actual and estimated currents of the third harmonic are impacted by the high frequency component of the PWM, especially when the amplitude of the current is low.

Even though the third harmonic of the FIPMSM (parameters provided in Table 2) only accounts for 10% of the first harmonic, the sliding mode observer allows a precise estimation of the back-EMF of the secondary fictitious machine in wide speed range. Thus, it can be concluded that the third harmonic component can be used, in the case of FIPMSM, to estimate the rotor position.

In the following section, actual rotor position, speed and torque will be compared to the estimated ones using the proposed sensorless control based on SMO.

Figure 5 and 6 show the actual and estimated values of the rotor position and speed. The simulations waveforms are obtained for the reference rotor speed and torque given by the Figure 2. In fact, the SMO allows an accurate estimation of the rotor speed in steady state and transient state as shown in Figure 5 (a). From Figure 5 (d), the estimated rotor position 𝜃 ̂𝑝 through the main fictitious machine converge to the actual one with high accuracy. From Figure 6 (b), the estimated rotor position 𝜃 ̂𝑠 through the second fictitious machine is also precise. Using the angles 𝜃 ̂𝑝 and 𝜃 ̂𝑠 for the sensorless control of the FIPMSM, it can be noticed that the torque converges to the actual one (Figure 5 (b)) obtained when the rotor position is provided by encoder.

The error between the real and estimated θp, through the main fictitious machine (first harmonic), is less than 1.5 degree as shown in Figure 5 (d). In addition, the error between the real and estimated θs, through the secondary fictitious machine (third harmonic), is less than 6 degrees as shown in Figure 6 (b). This error between the real and estimated position is evaluated in the medium speed range (100-1300 rpm). Therefore, it is concluded that, in the case of FIPMSM, the third harmonic can also be used to perform the sensorless control. However, at zero and low speed range (0-100 rpm), the rotor position estimation is not accurate as in Figures 5 (c)-(d), and Figures 6 (a)-(b). This is due to the low amplitude of the back-EMF at low speed range. In another hand, it should be noted that the pulses present in the position error (Figure 5 (d) and Figure 6 (b)), resulting to the compute of the estimation error, have not any impact on the control loop. This is due to the trigonometric functions used by the Park matrix transformation allowing the compute of currents and voltages components in (d-q) frame. From Figure 6 (c) and (d), the robustness of the proposed sensorless control based on SMO during speed reversal is verified. The rotor speed is estimated with accuracy and a good quality of torque under sensorless control mode is achieved.

Verifying Robustness of the proposed sensorless control

The robustness of the proposed observer, based on SMO, is required to perform an efficient sensorless control of the FIPMSM. In the above section 4.1, the robustness of the observer against speed and torque variation is verified. Thus, the robustness when the parameters of the FIPMSM (resistance, inductance, and flux of permanent magnets) change should be also verified. It is to highlight the efficiency of the proposed sensorless control approach.

The verification of the robustness against the machine parameters variation can be tested by using the defined benchmark given in [START_REF] Ghanes | Original benchmark for sensorless induction motor drives at low frequencies and validation of high gain observer[END_REF]. The reference rotor speed and the reference torque are shown in Figure 7. As defined in the benchmark [START_REF] Ghanes | Original benchmark for sensorless induction motor drives at low frequencies and validation of high gain observer[END_REF], the resistance as defined in the benchmark [START_REF] Ghanes | Original benchmark for sensorless induction motor drives at low frequencies and validation of high gain observer[END_REF], the resistance of the machine increases by 50% and decreases by 50%. The inductance increases by 20% and decreases by 20%, and the amplitude of the flux of permanent magnets increases also by 15% and decreases by 15%. This is to introduce a variation in the machine parameters. It is important to be noticed that the parameters variation is only at the machine (FIPMSM), and not on the control system and the SMO. It can be noticed that the machine parameters are affected mostly by the temperature, the saturation and the frequency [START_REF] Bogado | Sensitivity to electrical parameter variations of predictive current control in multiphase drives[END_REF]. The simulation results.

With variations of the resistance and the inductance of the FIPMSM are given respectively in Figures 8 and9.

In Figure 8, the resistance variation (due to the variation of the temperature inside the machine) has not an impact on the measured torque. However, the estimated speed when the resistance increases by 50% is not precise especially when the torque increases (Figure 8 (c)). The error estimation value between the actual and estimated rotor speed is 30 rpm when the torque is 10Nm. Thus, it can be concluded that the proposed SMO is robust and can guarantee the sensorless control under an important resistance variation. Simulation waveforms when the inductance of FIPMSM is changed: (a) the measured torque (using an encoder) and the one (using sensorless) when the inductance decreases by 20%, (b) the measured torque (using an encoder) and the one (using sensorless) when the inductance increases by 20%, (c) Actual and estimated rotor speed when the inductance decreases and increases by 20%

From Figure 9, the inductance variation (due to the saturation effects of magnetic circuit of the machine) presents no significant impacts on the measured torque. Nevertheless, the error estimation of the rotor speed is not precise when the torque is applied.

The error estimation is less than 20 rpm (Figure 9 (c)), when the inductance increases by 20% and decreases by 20%. So, based on the results in Figure 9, it can be concluded that the proposed SMO is robust and can guarantee the sensorless control under the inductance variation. Figure 10 shows results under the flux of permanent magnets variation (due to the variation of the magnets temperature). When the amplitude of flux decreases by 15%, the error between the measured torque using an encoder and the one using the proposed sensorless control is less than 6%.

And the error when the amplitude of flux increases by 15% is less than 4%. In Figure 10 (c), the estimated rotor speed is greatly affected in both cases, when the flux increases and decreases by 15%. This is due to the estimation process of the rotor speed, that depends directly on the amplitude of the flux as shown in [START_REF] Qiao | New sliding-mode observer for position sensorless control of permanent-magnet synchronous motor[END_REF]. Thus, it can be concluded that the robustness of the proposed SMO is affected under the flux variation of permanent magnets.

CONCLUSIONS

In this paper, an observer based on the sliding mode is designed for the sensorless control of a non-sinusoidal back-EMF of FIPMSM. The approach can be easily generalized to a n-phase machine whose electromotive force contains only the harmonics from one to (n-1)/2. On the contrary, the limit of the method is appearing for a five-phase (resp n-phase) machine if harmonic higher than the fifth (resp. the nth) are present. The simulation results show that the proposed SMO estimates the rotor position and speed in a wide speed range with high accuracy. Furthermore, in terms of robustness, the results show that the observer can guarantee the precise sensorless control under the resistance and inductance variation. However, the robustness is impacted under the flux of permanent magnets variation, especially in the rotor speed estimation process. Therefore, the proposed SMO with the consideration of the impact of third harmonic component (in FIPMSM) allows an accurate estimation of the rotor position, as the conventional SMO (with the consideration of only first harmonic). It can be noticed that when the third harmonic is important, the estimation of rotor position through the secondary fictitious machine is more precise. Thus, it can be concluded that the proposed observer can significantly improve the robustness and reliability of the control system by increasing the degrees of freedom for control. Lastly, the proposed SMO sensorless control in this paper does not take into account the low speed range (0-10% of the base speed). Indeed, the developed algorithm is applied to a 48V FIPMSM since 10% of speed leads to a very low back-EMF which cannot be estimated correctly.

13 

 13 𝑒 𝛽3 ][ 𝑖 𝛼3 𝑖 𝛽3] 𝑇 𝛺

Figure 1 .

 1 Figure 1. Block diagram of the proposed SMO for the sensorless control of FIPMSM with the consideration of third harmonic component

Figure 2 .

 2 Figure 2. (a) Reference rotor speed, (b) Load torque

Figure 3 .

 3 Figure 3. Simulation waveforms: (a) Actual and estimated fundamental current in α-axis, (c) Actual and estimated current of third harmonic in α-axis, (b) and (d) the error obtained from the estimation of fundamental and third harmonic of current

Figure 4 .

 4 Figure 4. Simulation waveforms: (a) Actual and estimated fundamental of back-EMF in α-axis, (c) Actual and estimated third harmonic of back-EMF in α-axis, (b) and (d) The error resulting from the estimation of fundamental and third harmonic of back-EMF

Figure 5 .Figure 6 .

 56 Figure 5. Simulation waveforms: (a) Actual and estimated rotor speed, (b) the measured torque (using an encoder) and the one (using sensorless) (c) sine actual and estimated rotor position of the main fictitious machine, (d) actual and estimated rotor position of the main fictitious machine

Figure 7 .

 7 Figure 7. Benchmark used to verify the robustness of the sensorless control: (a) reference rotor speed, (b) reference torque

Figure 8 .Figure 9 .

 89 Figure 8. Simulation waveforms when the resistance of FIPMSM is changed: (a) the measured torque (using an encoder) and the one (using sensorless) when the resistance decreases by 50%, (b) the measured torque (using an encoder) and the one (using sensorless) when the resistance increases by 50%, (c) Actual and estimated rotor speed when the resistance decreases and increases by 50%

Figure 10 .

 10 Figure 10. Simulation waveforms when the flux of the permanent magnets of FIPMSM is changed: (a) the measured torque (using an encoder) and the one (using sensorless) when the flux decreases by 15%, (b) the measured torque (using an encoder) and the one (using sensorless) when the flux increases by 15%, (c) Actual and estimated rotor speed when the flux decreases and increases by 15%

Table 1 .

 1 Fictitious machines and associated harmonics of FIPMSM[START_REF] Semail | Right harmonic spectrum for the back-electromotive force of an n-phase synchronous motor[END_REF]

	Fictitious machines	Associated harmonics
	Main machine	1, 9, 11, …5*k±1
	Secondary machine	3, 7, 13, …5*k±2
	Homopolar machine	5, 15, 25, …5*k

Table 2 .

 2 Parameters of the FIPMSM

	Parameters	Units	Values
	Rated voltage (Vdc)	V	48
	Rated power	kW	8
	Base speed	rpm	1300
	Speed-normalized amplitude of 1 st harmonic EMF	V/rad/s	0.1358
	Speed-normalized amplitude of 3 rd harmonic EMF	V/rad/s	0.01356
	ψ1, Magnetic flux of the 1 st harmonic	mWb	19.4
	ψ3, Magnetic flux of the 3 rd harmonic	mWb	0.675
	Stator resistance	mΩ	11
	Lp, Inductance of the main fictitious machine	µH	118
	Ls, Inductance of the second fictitious machine	µH	51.4
	Pole pairs		7

Table 3 .

 3 Parameters of the SMO for the sensorless control

	Parameters	k1	k2	l1	l2	a
	Values	250	25	500	1000	0.1

ACKNOWLEDGMENT

This work has been achieved within the framework of CE2I project. CE2I is co-financed by European Union with the financial support of European Regional Development Fund (ERDF), French State and the French Region of Hauts-de-France.