Hêmin Golpîra
email: golpira@wisc.edu

Azin Atarodi
email: a.atarodi@eng.uok.ac.ir

Shiva Amini
email: sh.amini@eng.uok.ac.ir

Fellow, IEEE Arturo Roman Messina
email: arturo.roman@cinvestav.mx

Hemin Golpîra
email: hemin.golpira@uok.ac.ir

Hêmin Golpîra

Hassan Bevrani
email: bevrani@uok.ac.ir

Bruno Francois
email: bruno.francois@centralelille.fr

Optimal Energy Storage System-Based Virtual Inertia Placement: A Frequency Stability Point of View 1 2

Keywords: Optimal placement, frequency nadir, virtual 28 inertia, energy storage systems, inertial response, rate of change 29 of frequency, transient stability, uncertainty analysis, sensitivity 30 analysis

In this paper, the problem of optimal placement of 5 virtual inertia is considered as a techno-economic problem from a 6 frequency stability point of view. First, a data driven-based equiv-7 alent model of battery energy storage systems, as seen from the 8 electrical system, is proposed. This experimentally validated model 9 takes advantage of the energy storage system special attributes to 10 contribute to inertial response enhancement, via the virtual inertia 11 concept. Then, a new framework is proposed, which considers the 12 battery storage system features, including annual costs, lifetime and 13 state of charge, into the optimal placement formulation to enhance 14 frequency response with a minimum storage capacity. Two well-15 known dynamical frequency criteria, the frequency nadir and the 16 rate of change of frequency, are utilized in the optimization formu-17 lation to determine minimum energy storage systems. Moreover, a 18 power angle-based stability index is also used to assess the effect 19 of virtual inertia on transient stability. Sensitivity and uncertainty 20 analyses are further conducted to assess the applicability of the 21 method. The efficiency of the proposed framework is demonstrated 22 on a linearized model of a three-area power system as well as two 23 nonlinear systems. Simulation results suggest that the proposed 24 method gives improved results in terms of stability measures and 25 less ESS capacity, when compared with other methods proposed in 26 the literature.

I NERTIAL response is defined as the power associated with changes in kinetic energy of synchronous generator rotors, in response to frequency changes, which is fed(taken) to(from) the grid [1]. Before reaction of traditional ancillary control loops, this energy, provided by system inertia, is the key factor to limit the power imbalance. System inertia has a major influence on frequency stability and the associated rate of change of frequency (RoCoF) and frequency nadir.

0885-8950 © 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information.

Inverter-based renewable sources are increasingly replacing synchronous generators, which in turn decrease overall system inertia. The ever-growing number of frequency incidents, in response to fluctuations of renewable power sources, accompanied with low level inertia feature jeopardize frequency stability

[2]- [4]. This motivates the need to develop advanced ancillary energy balancing services to control frequency changes.

Paramount among these, is Virtual Inertia emulation which mitigates undesired frequency dynamics. However, not only insufficient level of inertia, but also its heterogeneous distribution together with time-varying inertia profiles may render frequency dynamics faster. These facts, along with the need to economically keep system secure, make optimal placement of virtual inertia a key factor [2]. This paper addresses the problem of optimal placement of virtual inertia in power grids from a fundamentally new perspective.

Effects of energy storage systems (ESSs) on frequency regulation have been studied in recent research works such as [5]- [START_REF] Zhang | Assessment of the effectiveness of energy storage resources in the frequency regulation of a single-area power system[END_REF]. They deal with balancing of generation and load to maintain a constant system frequency and to keep tie-line power flows at scheduled values. These studies consider long term frequency response as well as steady state metrics, while neglect inertia requirements and primary frequency as the main metrics utilized for system resilience analysis. Some recent studies have investigated the effect of low system inertia on frequency stability [START_REF] Ahmadyar | A framework for assessing renewable integration limits with respect to frequency performance[END_REF]- [START_REF] Chu | Towards optimal system scheduling with synthetic inertia provision from wind turbines[END_REF]. In parallel with these efforts, recent research works [START_REF] Golpîra | Emulation of virtual inertia to accommodate higher penetration levels of distributed generation in power grids[END_REF]- [START_REF] Poolla | Placement and implementation of grid-forming and grid-following virtual inertia and fast frequency response[END_REF] demonstrate the way on which virtual inertia could be emulated in different ways, including appropriate control of wind turbines and energy storage systems (ESSs). In [START_REF] Borsche | Effects of rotational inertia on power system damping and frequency transients[END_REF], an optimization framework to deal with various aspects of inertia emulation and control including how inertia emulation impacts system stability, and determining the best places to add virtual inertia is proposed. Some questions about the heterogeneous inertial profiles and how the associated negative impacts are reduced by inertia emulation has been raised in [START_REF] Ulbig | Impact of low rotational inertia on power system stability and operation[END_REF]. Poolla et al.

[2] and Farmer et al. [START_REF] Farmer | Optimising power system frequency stability using virtual inertia from inverter-based renewable energy generation[END_REF] focus on H 2 performance metric to determine optimal placement of virtual inertia. Determination of the optimum size of battery to provide the primary frequency control is addressed in [START_REF] Fini | Determining optimal virtual inertia and frequency control parameters to preserve the frequency stability in islanded microgrids with high penetration of renewables[END_REF], [START_REF] Oudalov | Optimizing a battery energy storage system for primary frequency control[END_REF]. In particular, ESSs do not constantly participate in primary regulation, due to life-time concerns. Instead, in practice, the ESSs are controlled in a hierarchical fashion, dispatched by an optimal function to guarantee the State-of-Charge (SOC) and life-time constraints, which is missed in the so far researches. Some research works, such as [START_REF] Wogrin | Optimizing storage siting, sizing, and technology portfolios in transmission-constrained networks[END_REF], deal with optimal placement of virtual inertia in power system considering network structure. The main step towards optimal placement of virtual inertia in 160 a power grid is to analyze its effects on the frequency stability 161 behavior. This could be visualized by appropriate modelling 162 of the virtual inertia. Generally, providing virtual inertia and 163 thereby contributing to the overall equivalent grid inertia could 164 be achieved by using the virtual synchronous generator (VSG) 165 concept. VSG relies on similar power-balance-based synchro-166 nization mechanism as SGs to realize such functionality [START_REF] Mo | Evaluation of virtual synchronous machines with dynamic or quasi-stationary machine models[END_REF], 167 which would be investigated for modeling purpose in this paper. 168 With synchronous generation being displaced with modern 169 power electronic based generation such a solar and wind gen-170 eration, the analysis of inertial response of variable resource 171 generation becomes of fundamental importance. Energy storage 172 systems can, in principle, provide most of its stored energy 173 to support frequency in an interconnected power system and 174 hence a set of large battery ESSs could play the same role as 175 SGs in the inertial response horizon. Conceptually, the charg-176 ing/discharging process of ESS can be interpreted as the initial 177 compensation of a disturbance by the stored kinetic energy of 178 a SG rotating mass. This means that, while the ESSs would be 179 triggered by a command signal, from frequency measurements, 180 for charging/discharging, a SG model could be used to represent 181 the associated dynamics.

182

The classical model of a SG is a second-order model of the 183 form

184 δs = ω s (t) -ω 0 H ωs = T m (t) -T e (t) (1)
where δ s , ω s , and ω 0 are the mechanical rotor angle, the mechan-185 ical rotor angular speed and the initial angular speed, respec-186 tively; H, T m (t) and T e (t) are the inertia constant, mechanical 187 input torque and electrical output torque, respectively [START_REF] Ajala | A library of second-order models for synchronous machines[END_REF].

Taking the slow electromechanical behavior of the battery ESS into account, the associated dynamics could be represented by (1). The problem of interest, however, is to calculate the equivalent inertia constant and mechanical input torque. To deal with this possibility, a data-driven approach in which uncertain behavior of the ESS is inherently considered to provide complementary information for the swing equation model of (1), is proposed.

B. Data Driven-Based Equivalent Model:

To introduce the proposed equivalent model, assume that the injected power of the ESS to the host grid is a discrete-time signal P(n) of length L. Let the time-varying signal P(n) be decomposed into K sinusoidal components in noise, as [START_REF] Bollen | Signal Processing of Power Quality Disturbances[END_REF], [START_REF] Manolakis | Statistical and Adaptive Signal Processing: Spectral Estimation, Signal Modeling, Adaptive Filtering, and Array Processing[END_REF]:

P (n) = K k=1 a k cos (nω k + Φ k ) + w (n) (2)
where, a k , Φ k , ω k and w(n) are the magnitude and the initial phase angle, harmonic frequency in radius and additive white noise, respectively. In the model, a k and ω k are assumed to be deterministic and unknown, and Φ k is unknown and assumed to be random and uniformly distributed in [-π, π]. Alternatively, the model (2) can be expressed in the form of noisy complex exponentials as [START_REF] Farmer | Optimising power system frequency stability using virtual inertia from inverter-based renewable energy generation[END_REF]

P (n) = K k=1 A k e (jnω k ) + w (n) (3)
where A k = |A k |e jΦ k is the complex magnitude of the kthharmonic(noise) signal component. As the MUSIC algorithm is a noise subspace-based method, it is a good tool to deal with experimental noisy measured signals. Using this framework, the dimensional space is divided into signal and noise components, which is of high importance to accurately calculate H and T m (t)

in (1).

The MUSIC method employs a harmonic model and estimates the frequencies and powers of the harmonics in the signal.

Application of the MUSIC method to the data sequence, P(n),

gives:

P M U SIC e jω = 1 M i=K+1 |e U s i | 2 (4)
where the s i are the eigenvectors associated with the noise subspace that are orthogonal to the signal eigenvector e =

[1 e jω ej 2ω e j(M -1)ω ] T , and e U denotes the complex-conjugate transpose; M is the dimension of space spanned by P(n). It is worth emphasizing that P M U SIC (e jω ) in (4) does not relate to any real power spectrum; rather, the only purpose of this pseudospectrum is to generate peaks whose frequencies correspond to those of the dominant frequency components. This feature makes the MUSIC approach interesting to develop equivalent model based on dominant modes.

For a given signal of interest and according to (2)-(4), eigenvalues would be calculated. By knowing the eigenvalues and because the impulse response is the inverse Laplace transform of eigenvalues, one could represent a signal of interest with a pre-defined model of (1). For this purpose, suppose the im-234 pulse response of system is I(t); for the input signal x(t), i.e. 235 x(t) = T m (t) in (1), one could write y(t), i.e. ω, as:

236 y (t) = x (t) ∞ t=0 I (t) e jωt dt (5.a)
where, 237

I(t) = y(t) * P S (5.b)
and PS in (5.b) is the pseudo-spectrum of the signal. Equation, 238 (5.b) reveals that I(t) results from convolution of y(t) and PS.

239

In the modelling procedure and by measuring the output 240 response of the system y(t) and by knowing I(t), the problem 241 of interest is to calculate x(t) in (5). By calculating T m (t), the 242 ESS could be replaced by the synchronous generator model of 243 (1). Using this approach, o in (5.a) is defined as the dominant 244 frequency components of the pseudo-spectrum in (4). The input to the control algorithm is the frequency at the con-250 nection point of the inverter f cp , and P in represents the grid 251 injected power.

252

C. Model Validation:

The efficiency of the proposed method 253 is illustrated using the existing battery ESS in the University of 254 Kurdistan Micro-Grid (UOK-MG). Fig. 2 shows a three-phase 255 diagram representation of the UOK-MG. Details of the physical 256 UOK-MG are given in [START_REF] Golpîra | Bulk power system frequency stability assessment in presence of microgrids[END_REF].

257

Figure 3 shows the battery ESS and the main grid power 258 variation behavior, i.e. P in in Fig. 1, recorded for the UOK-MG. 259 As shown in Fig. 3, Event 1 triggers the charging process of the 260 ESS in response to deviations from the minimum SOC.

261

The main grid power deviation during the charging process in 262 Fig. 3 is utilized to calculate the Pseudo-spectrum (see Fig. 4) 263 which, in turn, is used to estimate the dominant frequency 264 components in (4). Setting the frequency deviation of the ESS in Fig. 5 as y(t) 266 and the Pseudo-spectrum of Fig. 4 as PS for (5) gives:

267 δs = ω s -ω 0 0.53 ωs = 1 -e -0.38t C -T e (t) (6)
where, C is a constant value and would be set to fit the DC term in 268 (4). In interpreting this model, note that constant 0.38 in (6) rep-

269

resents the dominant frequency of the Pseudo-spectrum in Fig. 4 270 which would be fed to (5) to calculate the 0.53 s inertia constant in (6). Figure 5 demonstrates the effectiveness of the equivalent model (6) to approximate the inertial response behavior of the ESS. To exactly mimic frequency behavior of ESS using (6), the oscillatory behavior of the adopted model would be removed using a 20-samples rolling-averaging window. This approach averages the long-term oscillations, and hence, mitigates the oscillatory behavior beyond the inertial response horizon.

Results in Fig. 5 show that the dynamic behavior of ESS, especially in the inertial response horizon can be approximated by a simple SG model. It should be emphasized that while a conventional SG is slower and less flexible compared to ESS, developed dynamical equivalent model of ESSs in the grid connected mode is not only affected by the fast-inherent features of ESSs but also is significantly influenced by the external zone, i.e. the host grid, features.

III. OPTIMIZATION PROBLEM

This section formulates the ESS placement as an optimal techno-economic problem.

A. Costs of Energy Storages and Technologies

There are two main approaches for assessing the cost of storage technologies: 1) Total Capital Cost (TCC), and 2) Life Cycle Cost (LCC) [START_REF] Zakeri | Electrical energy storage systems: A comparative life cycle cost analysis[END_REF]. For the sake of generality of the method, no specific ESS technologies are considered in this section.

In the TCC approach, all terms associated with the purchase, installation, and delivery of the ESS units, including Power Conversion System (PCS) costs (C P CS ), costs of ESS (C stor ), and cost of balance of plant (C BOP ), are considered as:

C cap = C P CS + C BOP + C stor × t ch ($/kW ) (7)
where t ch is the charging /discharging time. The balance of the ESS, known as the BOP, includes site wiring, interconnecting transformers, and other additional ancillary equipment and is measured on a $/kW basis [START_REF] Mongird | Energy storage technology and cost characterization report[END_REF].

However, the LCC is a more common metric to evaluate and compare different ESS technologies. The annualized LCC is formulated according to (8) which considers operation and maintenance costs (C o&M,a ), replacement cost (C R,a ) and annualized TCC.

C LCC,a = C cap,a + C O&M,a + C R,a ($/kW -yr) (8)
in which:

CRF = i(1 + i) T (1 + i) T -1 (9) C cap,a = T CC × CRF ($/ (kW -yr)) (10)
C O&M,a = C F OM,a +C V OM,a × n cycle × t ch ($/kW -yr) (11) C R,a = CRF × r k=1 (1 + i) -kt × C R × t ch η sys ($/kW -yr) (12)
where CRF, i, T, r, t, and η sys are the capital recovery factor, interest rate and the life time, the number of substitutions in lifetime, the replacement period and the overall efficiency, respectively; C F OM,a and C V OM,a define the fixed and variable operation and maintenance costs. Subscript "a" stands for "annualized" costs.

B. Formulation of Objective Function and Constraints

Equation ( 8) specifies the annual cost per kilowatt of the installed ESS in compliance with the lifetime. However, for optimal placement of virtual inertia, it is necessary to rewrite the cost function in [START_REF] Ahmadyar | A framework for assessing renewable integration limits with respect to frequency performance[END_REF] according to the amount of virtual inertia.

The synchronous inertia constant (H) is defined as the ratio of stored kinetic energy to the rated apparent power of the system, as:

H = 0.5J V I ω 2 S base (13)
where J V I , ω, and S base are the moment of inertia, angular velocity and rated apparent power, respectively. Since the stored energy in ESSs is usually expressed in volt ampere hour (V Ah ESS ), it is needed to express the associated value in term of Joule. Considering unity power factor, one could re-write (13) as: calculate inertia constant as [START_REF]Fast Freq. Response Concepts Bulk Power Syst. Rel. Needs[END_REF]:

KW ESS = KV A ESS = h ESS S base 3600 s → h ESS = 3600 V Asec ESS S base (
H Conv. = ∆P L 2RoCoF 500ms = ∆P L 2 f (0.5)-f (0) 0.5 (14.b)
A comparison with the frequency responses of Fig. 5 shows an error of 3.17% which in turn justifies the proposed formulation in (14.a).

By substituting (14.a) into (8), one could write the optimization problem as: The SOC can be calculated as follows [START_REF] Teruo | State of charge calculation device and state of charge calculation method[END_REF][START_REF] Tang | Frequency regulation mechanism of energy storage system for the power grid[END_REF] 343

minimize h i F (h i ) = n ESS i=1 C LCC,a i × h i S
SOC (∆t) = SOC (0) - ∫ ∆t 0 ζp (t) dt E ESS,rated (16.

1) Rate of Change of Frequency:

The RoCoF is a meaningful 360 criterion to show ability of a system in the face of a sudden 361 power imbalance. Greater RoCoF means that the less time is 362 available for system operator to arrest frequency decline [START_REF] Eto | Use of frequency response metrics to assess the planning and operating requirements for reliable integration of variable renewable generation[END_REF]. 363 Time interval of 100 milliseconds to 2 seconds are defined to 364 measure the RoCoF [START_REF] Eto | Use of frequency response metrics to assess the planning and operating requirements for reliable integration of variable renewable generation[END_REF], [START_REF]Frequency Stability Evaluation Criteria for the Synchronous Zone of Continental Europe -Requirements and impacting factors[END_REF]. ENTSO standard [START_REF]Frequency Stability Evaluation Criteria for the Synchronous Zone of Continental Europe -Requirements and impacting factors[END_REF] explains 365 that RoCoF is allowed to get a value between 0.5 to 1 Hz/sec.

366

In order to represent dynamical frequency indices based on 367 lower bounds of inequality constraints of (15.b), the RoCoF 368 would be defined based on the classical swing equation of (1) 369 as [START_REF] Golpîra | Bulk power system frequency stability assessment in presence of microgrids[END_REF]:

370 2H d∆f (t) dt = ∆P m (t) -∆P L (t) -∆P tie (t) (17)
Where ∆P m (t), ∆P L (t), and ∆P tie (t) represent mechanical 371 power, electrical power and tie line power changes, respectively. 372 Considering the definition of RoCoF, one could write:

373 RoCoF = ∆P m (t) -∆P L (t) -∆P tie (t) 2H (18)
The Taylor series expansion of (18) about the independent 374 variables of H, ∆P m , ∆P L , and ∆P tie gives

375 ∆RoCoF i = ∂RoCoF i ∂∆P mi ∆∆P mi + ∂RoCoF i ∂∆P Li ∆∆P Li + ∂RoCoF i ∂∆P tiei ∆∆P tiei + ∂RoCoF i ∂∆H i ∆H i = 1 2H i ∆∆P mi + -1 2H i ∆∆P Li + -1 2H i ∆∆P tiei + -(∆P mi -∆P Li -∆P tiei ) 2H i 2 ∆H i (19)
Due to the slow inherent of dynamics of interest, except for the 376 last term of ( 19), other terms could be neglected. Accordingly, 377 and by replacing [START_REF] Hammad | On effective virtual inertia of storage-based distributed control for transient stability[END_REF], one could re-write [START_REF] Farmer | Optimising power system frequency stability using virtual inertia from inverter-based renewable energy generation[END_REF] as:

378 RoCoF i - ∆H i H i = ∆RoCoF i (20)
Considering maximum allowable RoCoF, i.e. RoCoF i,max , 379 the minimum inertia which guarantees the RoCoF to be within 380 the permitted range is calculated as:

381 h ′ i,min = ∆H i,min = H i - ∆RoCoF i,max RoCoF i ∆RoCoF i,max =RoCoF i,max -RoCoF i ------------------------→ h ′ i,min = H i - RoCoF i,max -RoCoF i RoCoF i (21)
where

393

Taking the time dependency of the governor response into 394 account, one can write the frequency nadir as [START_REF] Golpîra | Emulation of virtual inertia to accommodate higher penetration levels of distributed generation in power grids[END_REF]:

395 ∆f nadir = (∆P L + ∆P tie ) 2 T d 4HR ( 22
)
where R is the extra power received through the governor and 396 T d is the response time of the governor. In deriving [START_REF] Borsche | Effects of rotational inertia on power system damping and frequency transients[END_REF]

+ ∂∆f nadir,i ∂∆H i ∆H i = (∆P Li + ∆P tiei ) T di 2H i R i ∆∆P Li + (∆P Li + ∆P tiei ) T di 2H i R i ∆∆P tiei + -(∆P Li + ∆P tiei ) 2 T di 4H 2 i R i ∆H i ( 23
)
Following the same procedure as that in [START_REF] Poolla | Placement and implementation of grid-forming and grid-following virtual inertia and fast frequency response[END_REF], one could rewrite [START_REF] Fini | Determining optimal virtual inertia and frequency control parameters to preserve the frequency stability in islanded microgrids with high penetration of renewables[END_REF] in the form

∆f nadir,i - ∆H i H i = ∆∆f nadir, i = ∆(f nadir, i -f 0 ) (24)
The minimum inertia, i.e. h ′′ i,min , which guarantees frequency nadir to be in the permitted range is calculated by:

h ′′ i,min = ∆H i,min = H i - ∆∆f nadir i, max ∆f nadir i = H i - f nadir max -∆f nadir i ∆f nadir i ( 25
)
where, h ′′ i,min represents the minimum, in compliance with frequency nadir, required inertia which should be emulated by the battery ESS in area i. In order to simultaneously satisfy both, frequency nadir and RoCoF standards, the lower bound for virtual inertia in the optimization problem and for each area are selected as the maximum value of ( 21) and ( 25), namely:

h i,min = max h ′ i,min , h ′′ i,min (26)
Moreover, the overall system inertia has a direct impact on the frequency indices. This means that some considerations should be made regarding overall system inertia and, consequently, ( 21) and ( 25) would be completed by adding a new equality constraint. For this purpose, the frequency of the overall Center of Inertia (COI), which should satisfy strict frequency standards, would be employed to determine the overall amount of inertia in the system, as [5]:

H COI = Q = H ∆f COI f COI (61.5) (27)
where, ∆f COI , and f COI represent the frequency deviation and frequency of the system, without ESS, after the fault, respectively. Formally, equation [START_REF] Ajala | A library of second-order models for synchronous machines[END_REF] gives the required amount of inertia constant which guarantees acceptable frequency dynamics of the COI. Of note that Q would be realized by adding the emulated inertia of ESSs to the conventional SGs inertia. Accordingly, the optimization problem ( 15) can be re-written as:

minimize h i F (h i ) = n ESS i=1 C LCC,a i h i S base 3600 sec (28.a)
st:

H COI = Q (28.b) h i,min ≤ h i ≤ h i,max (28.c) SOC min ≤ SOC i ≤ SOC max (28.d)
where the dynamic inequality constraints (15.b) and (15.c) are re-formulated as the algebraic inequality constraint (28.c) in terms of the inertia constant. This dramatically increases the simplicity and speed of the calculations.

IV. SOLUTION ALGORITHM

In this section, a framework that incorporates the model (6) and the objective function [START_REF] Bollen | Signal Processing of Power Quality Disturbances[END_REF] to optimally place battery ESS in the system is proposed. The algorithm consists of 4 main steps:

Step 1: Define PV buses of the system as candidates for battery

ESSs placement;

Step 2: Assume the proposed model ( 6), with unknown inertia constant, in each PV bus;

Step 3: Define the emulated inertia constant, i.e. h, as a decision variable of the objective function [START_REF] Bollen | Signal Processing of Power Quality Disturbances[END_REF]. Set C in (6) as the product of the ratio of h to the inertia constant of the installed SG, in the associated bus, and the mechanical power of the SG;

Step 4: Apply a Genetic Algorithm (GA) to solve the minimization problem [START_REF] Bollen | Signal Processing of Power Quality Disturbances[END_REF]. Indeed, the GA determines the amount of virtual inertia in each PV bus of the system.

V. SIMULATION AND RESULTS

Three test systems have been used for evaluating the proposed formulation in this paper: a) a linearized model of a three-area power system, b) a nonlinear yet simple two-area, four-machine test system, and c) a large scale16-machine, five-area 68-bus test model of the New York/New England system.

A. Linear System

As a first motivating example, a linearized model of a threearea power system is used to assess the efficiency of the proposed formulation. The block diagram of each area is shown in Fig. 6. Firstly, a 0.2 per unit load disturbance is applied at areas 1 and 3. As a first scenario, the required virtual inertia is calculated only based on (28.b) and arbitrary realized, through the proposed model in (6), in area 1. For simulation purposes, it is assumed that the ESS would be triggered upon occurrence of the fault.

Comparison of frequency dynamics for the system with and without virtual inertia reveals that while areas 1 and 2 frequency nadirs improved by means of virtual inertia, area 3 shows an undesired behavior. This in turn numerically justifies the need for optimal inertia placement. Within this framework, (28) may be written as: where, for instance, the minimum inertia of area 1 in (29.c) is 473 calculated based on ( 21), [START_REF] Wogrin | Optimizing storage siting, sizing, and technology portfolios in transmission-constrained networks[END_REF], and (26) as:

minimize h i F (h i ) = n ESS i=1 C LCC,a i h i S base 3600
474 h ′ 1,min = 0.08335 × - 1 + 1.1870 -1.1870 = 0.0129 h ′′ 1,min = 0 h 1,min = max {0, 0.0129} = 0.0129 (30)
Note that some other unspecified parameters of (29) would 475 be set as described in Appendix A. The results obtained, using 476 a simple GA with 0.05 and 0.8 mutation and crossover coef-477 ficients, respectively, from optimization of (29) are shown in 478 Table I. In order to further assess efficiency of the proposed 479 formulation, Table I compares the results with those of ob-480 tained from [2] and [START_REF] Borsche | On placement of synthetic inertia with explicit 785 time-domain constraints[END_REF]. Comparison results justify the fact 481 that dynamic behavior of ESS could significantly affect optimal 482 placement problem which is missed in previous research.

Figure 7 compares the frequency behavior and RoCoF of 484 generators 1, 2 through 3 for three cases of interest: a) with 485 virtual inertia and according to the proposed formulation, b) with 486 virtual inertia and according to [START_REF] Borsche | On placement of synthetic inertia with explicit 785 time-domain constraints[END_REF], and c) with virtual inertia 487 and according to [2]. It can be seen that while frequency traces 488 Also of interest, modal analysis of the results, as explained in Table II, shows the efficiency of the proposed formulation in comparison with that of [2] and [START_REF] Borsche | On placement of synthetic inertia with explicit 785 time-domain constraints[END_REF].

Results show that a lower emulated virtual inertia in the proposed formulation not only decreases the cost function but also causes better performance in terms of enhanced damping.

With virtual inertia, the results from [2] and [START_REF] Borsche | On placement of synthetic inertia with explicit 785 time-domain constraints[END_REF] also seem to perform within the constraints and are almost the same as the results from the proposed method. This could be justified through the fact that set of PV buses for small systems includes a few members to be considered as candidates of ESS installation.

Therefore, different algorithms may differ a bit from capacity point of view rather than location which in turn causes negligible difference between the results. To further assess simultaneous effects of virtual inertia on frequency and transient stabilities, two non-linear system are used in what follows.

B. Two-Area Power System

The two-area power system, shown in Fig. 8, is considered to further demonstrate the efficiency of the proposed formulation.

Modeling considerations are essentially those described in [START_REF] Rogers | Power System Oscillations[END_REF]; all the generating units are modeled with 6 th order synchronous machine models with excitation systems.

The disturbance of interest is the outage of generator G 4 , at the first second of simulation. The lower bounds of virtual inertia are calculated according to [START_REF] Manolakis | Statistical and Adaptive Signal Processing: Spectral Estimation, Signal Modeling, Adaptive Filtering, and Array Processing[END_REF] and [START_REF] Golpîra | Bulk power system frequency stability assessment in presence of microgrids[END_REF] as: III. The 518 results also compared with those of optimum virtual inertia of 519 [2] and [START_REF] Borsche | On placement of synthetic inertia with explicit 785 time-domain constraints[END_REF].

minimize h i F (h i ) = n ESS i=1 C LCC

520

Further, Table IV compares frequency stability indices for 521 different approaches. The results demonstrate high efficiency of 522 the proposed method to optimally allocate virtual inertia in the 523 system.

524

While the frequency nadir of the generators are within the 525 permissible range for the system without virtual inertia, the 526 RoCoFs for the generators 1 and 2 exceed the standard value. 527 Optimal placement of virtual inertia returns generators with 528 undesired frequency dynamics to the normal region. Efficiency 529 of the proposed method is further assessed through time domain 530 simulations of Fig. 9.

531

Also of interest, effects of virtual inertia on transient stability 532 is assessed using a simple power angle-based stability index 533 η [START_REF] Hammad | On effective virtual inertia of storage-based distributed control for transient stability[END_REF]:

Q3 534 η = 360 -δ max 360 + δ max ( 33
)
where δ max is the maximum angle separation between any two 535 generators in the system [START_REF] Inc | Transient security assessment tool user manual[END_REF]. Generally stated, the use of simple 536 metrics such as [START_REF]Fast Freq. Response Concepts Bulk Power Syst. Rel. Needs[END_REF]

C. New York/New England System

The New England test system is used to further illustrate the efficiency of the proposed algorithm for large scale power systems. A single-line diagram of the system, showing major coherent areas and their interconnections, is shown in Fig. 10.

Five different contingency scenarios, including tripping of major generating units and load shedding, are considered. Results in Table VI compare the proposed algorithm results with those of [2] and [START_REF] Borsche | On placement of synthetic inertia with explicit 785 time-domain constraints[END_REF]. The results show high efficiency of the proposed method to improve frequency dynamics with minimum cost. Also of interest, Fig. 11 shows the allocation of virtual inertia among the PV buses of the system. Moreover, Fig. 12 compares the frequency behavior of the system, in response to the outage of generator n1 (scenario 1), for the proposed method, [2] and [START_REF] Borsche | On placement of synthetic inertia with explicit 785 time-domain constraints[END_REF]. It should be noted that while there are negligible deviations between the traces in Fig. 12, significant differences between the cost functions justify the efficiency of the proposed method.

Results show that the optimization problem works better and more effectively for larger areas. This can be understood by noting that set of PV buses for large areas includes many members to be considered as candidates for ESS placement.

As a result, there are a lot of possibilities for both the size and location. Moreover, the efficiency of the proposed formulation to en-574 hance transient stability is shown in Table VII, using [START_REF]Fast Freq. Response Concepts Bulk Power Syst. Rel. Needs[END_REF].

575

Table VII shows that appropriate placement of virtual inertia 576 in the system, considering dynamical behavior of ESS, could 577 also improve transient stability. This could be justified through 578 the fact that the so far researches, e.g. [2], rely on quasi-steady 579 state phasors for voltages and currents in transient stability 580 assessment. In other words, they consider constant nominal 581 value of frequency in defining system impedances which is far 582 from realistic for system with penetration of inverter based ESSs. 583 This point is successfully addressed in the proposed formulation 584 by explicitly representing the dynamic behavior of ESSs in the 585 problem formulation.

F (h i ) + ∆F (h i ) = n ESS i=1 C LCC,a i h i S base 3600 sec + ∂F (h i ) ∂∆P L ∆∆P L + ∂F (h i ) ∂∆P m ∆∆P m + ∂F (h i ) ∂C LCC, a i ∆C LCC, a i (34)
Then, according to [START_REF] Farmer | Optimising power system frequency stability using virtual inertia from inverter-based renewable energy generation[END_REF]

The usefulness of ( 36) is now assessed for the New-England 597 system. For this purpose, contingency 1, i.e. the outage of 598 generator 1 in area 1, is considered as the base case. Cost function 599 for outage of generator 7 in area 1, i.e. scenario 2, and generator 600 11 in area 2, i.e. scenario 3, are calculated using [START_REF] Eto | Use of frequency response metrics to assess the planning and operating requirements for reliable integration of variable renewable generation[END_REF]. compares the results which clearly justify effectiveness of the proposed sensitivity analysis.

For uncertainty analysis, the equality constraint of (28.b) is represented in the objective function (28.a), as: minimize

h i F (h i ) = n ESS i=1
C LCC,a i h i S base 3600 sec

+ β(H COI -Q) ( 37
)
where b is arbitrary chosen high with the aim of enforcing the results to follow the equality constraint of (28.b). Considering parametric uncertainty for inertia constant i.e. H COI , one could write (37) as:

minimize h i F (h i ) = n ESS i=1
C LCC,a i h i S base 3600 sec

+ β(H COI + γ -Q) → minimize h i F (h i ) = n ESS i=1
C LCC,a i h i S base 3600 sec

+ β(H COI -Q) + βγ (38)
where g expressed in percentage of H COI . To deal with uncertainty analysis, a simple interval approach is utilized. It assumes that the uncertain parameters take value in a specified interval. It could be reinterpreted as the probabilistic modeling with a uniform probability density function (PDF). In this method, the upper and lower bounds for the uncertain inertia parameter are defined. The aim is to find the lower and upper bounds of objective function [START_REF] Aien | A comprehensive 791 review on uncertainty modeling techniques in power system studies[END_REF].

Using the proposed framework, assume that the maximum uncertainty of inertia constant is considered to be 5p. This means that the interval of interest can be defined as: γ = [H COI -0.05(H COI ), H COI + 0.05(H COI )] (39) which in turn causes F (h i ) as:

F (h i ) = [3.8298, 5.5112] (40)
with a uniform PDF.

VII. CONCLUSION

While rotational inertia stabilizes the frequency of power grids against small and large disturbances, it leads for oscillations between generators. This paper provides a framework for optimal placement of virtual inertia in low inertia power systems which in turn improves host grid frequency stability. In this way, the ESSs are used to emulate virtual inertia placement. On the other hand, the proposed method in this paper tackles dynamical behavior of the ESSs into problem formulation and thus causes less, in comparison with the literature, virtual inertia 632 to be implemented in the system. This in turn causes better rotor 633 angle stability features. Of note that the proposed algorithm can 634 be a valuable tool in generation expansion planning of power 635 system and inertia deployment.

636

A data driven-based approach to represent the ESS dynamics 637 using conventional synchronous generator is proposed. This in 638 turn causes the gathered data, from the field setup, provides

Fig. 1 .

 1 Fig. 1. Block diagram representation of the proposed modeling process: (a) virtual inertia emulation mechanism, and (b) the proposed equivalent model; P in , y(t) and f cp represent the grid injected power, the frequency deviation of the ESS, and the frequency at the point of connection of the ESS, respectively.

245Figure 1

 1 Figure 1 gives a schematic illustration of this model. In this 246 plot, Fig. 1(a) illustrates the process of virtual inertia emulation 247 using the battery ESS, while Fig. 1(b) describes a simplified 248 block diagram representation of the proposed equivalent model. 249The input to the control algorithm is the frequency at the con-250 nection point of the inverter f cp , and P in represents the grid 251 injected power.

Fig. 2 .

 2 Fig. 2. Three-phase schematic representation of the UOK-MG.

Fig. 3 .

 3 Fig. 3. The ESS and grid experimental dynamic responses.

Fig. 4 .

 4 Fig. 4. The Pseudospectrum estimation via MUSIC.

Fig. 5 .

 5 Fig. 5. Comparison of frequency response computed from the experiment and the equivalent frequency response models.

  14.a) This equation gives the average hourly power that can be injected/absorbed to/from the grid by the ESS. To validate such representation, results obtained from (14.a), for frequency response of Fig. 5, are compared with those of well-established methods of calculating inertia. Using the classical swing equation during 500 ms after fault occurrence, one can conventionally

Fig. 6 .

 6 Fig. 6. Block diagram representation of control area i. β i , R i , T gi , T ti , and D i are frequency bias, droop characteristic, governor time constant, turbine time constant, and damping property, respectively.

  sec (29.a) st: H COI = 0.053 (29.b)

Fig. 7 .

 7 Fig. 7. (a) Frequency behaviors and (b) RoCoF of generators 1, through 3 of three area power system for three cases: proposed formulation, Ref. [44] and Ref. [2].

Fig. 8 .

 8 Fig. 8. Single-line diagram of two area system.

Fig. 9 .

 9 Fig. 9. Frequency response of generators 1, 2 and 3 of the two area power system for three cases: proposed formulation, Ref. [44] and Ref. [2].

Fig. 10 .

 10 Fig. 10. Single line diagram of the 68-bus system showing coherent areas and their interconnections.

Fig. 12 .

 12 Fig. 12. Frequency responses of New England system for different approaches.

  An attempt is made to systematically represent the ESS dy-141 namical behavior by a second-order SG model that extends 142 previous work by the authors [15] on the emulation of vir-143 tual inertia. Within this framework, a new signal processing 144 technique, called MUltiple Signal Classification (MUSIC) 145 algorithm is developed to derive the model. As formulating differential equations in the optimization 147 problem is a hard task, a mathematical framework is pro-148 posed to represent dynamical metrics by algebraic inequal-149 ity constraints. In particular, the present paper describes the 150 complex optimization problem of virtual inertia placement 151 by a new simple formulation. Since, the ESSs sizing is a difficult problem in practice, and 153 considering the energy capacity, price, life-time and SOC 154 may not be addressed by the proposed synchronous-based 155 model; an attempt is made to tie the virtual inertia to the 156 ESS features.

		157
	II. VIRTUAL INERTIA EQUIVALENT MODEL	158
	A. General Concept	159
	It employs DC-power
	flow to tackle network structure into formulation.
	A main limitation of these approaches is the reliance on steady
	state considerations. It is well known, however, that dynamical
	frequency indices, such as RoCoF and frequency nadir, are
	important parameters in assessing frequency stability and the
	activation of protection schemes. This paper re-formulates the

virtual inertia placement problem in term of dynamical metrics. Using this approach, the effects of virtual inertia on the frequency indices can be assessed using a dynamic equivalent model obtained by mapping the electromechanical behavior of ESSs onto a second-order synchronous generator (SG) model.

The main contributions of this paper can be summarized as 139 follows: 140 r 146 r 152 r

  is battery power which gets negative values for the 345 charging procedure and positive values for the discharging pe-346 riod; E ESS,rated , ∆t, ζ c , and ζ d are the nominal energy capacity, 347 charge/discharge time, charging and discharging efficiencies of 348 the battery, respectively. .c) explain that the optimization 351 problem (15) enforces the RoCoF and frequency nadir in all 352 areas to be less than standard values. These terms make the 353 optimization problem difficult to deal with as it depends on 354 dynamical indices. Generally, it is common to specify the 355 lower/upper bounds based on different criteria, including capac-356 ity of equipment or budget. Therefore, the problem of interest 357 here is to re-write the upper and lower bounds of (15.b) and 358 (15.c) in terms of the emulated inertia, i.e. h.

			a)
	where,		344
	ζ =	ζ c P (t) < 0 1 ζ d P (t) > 0	(16.b)
	and p(t) 359

349 C. Determining the Bounds of Virtual Inertia 350 Constraints (15.b) and (15

TABLE I OPTIMIZATION

 I RESULTS IN THREE AREA SYSTEM

a.PM: Proposed Method

TABLE II MODAL

 II ANALYSIS OF THE SYSTEMS WITH VI

TABLE III OPTIMIZATION

 III RESULTS IN TWO AREA SYSTEM

TABLE IV FREQUENCY

 IV INDICATORS OF TWO AREA SYSTEM BEFORE AND AFTER APPLICATION OF OPTIMAL INERTIA VALUES

  may work for some systems but fail for some 537 others. During severe faults, most ESSs, if remain connected and 538 continue to inject active power which is not the case in many 539 regions, get saturated and cannot follow the frequency properly. 540 Using[START_REF]Fast Freq. Response Concepts Bulk Power Syst. Rel. Needs[END_REF] in this paper relies on a conservative assumption that 541

TABLE VI OPTIMIZATION

 VI RESULTS IN NEW YORK/NEW ENGLAND SYSTEM Fig. 11. Virtual inertia allocation for New-England system.

  and[START_REF] Fini | Determining optimal virtual inertia and frequency control parameters to preserve the frequency stability in islanded microgrids with high penetration of renewables[END_REF], one could write the sensi-

	594
	595	tivity matrix (26) as:
			∂h ′ i,min ∂h ′′ i,min	=	∂RoCoF i ∂∆P Li ∂∆f nadir,i ∂∆P Li	∂RoCoF i ∂∆P mi 0	∆∆P Li ∆∆P mi	(35)
	596	Considering (35) in (34) gives:
		∆F (h i ) =	∂F (h i ) ∂h i	∂h i ∂∆P L	∆∆P L +	∂F (h i ) ∂h i	∂h i ∂∆P m	∆∆P m
		=	n ESS i=1	C LCC,a i	S base 3600 sec	× max	∂RoCoF i ∂∆P Li	∆∆P L
			+	∂RoCoF i ∂∆P mi	∆∆P m ,	∂∆F nadir,i ∂∆P Li	∆∆P Li
		+	n ESS i=1		h i S base 3600 sec	∆C LCC, a i

Table VIII

 VIII

	TABLE VIII
	SENSITIVITY ANALYSIS
	601


Henri-Claude Nataf
email: henri-claude.nataf@univ-grenoble-alpes.fr

H.-C Nataf

•••-• Tidally synchronized solar dynamo: a rebuttal

Keywords: Solar Cycle, Models; Solar Cycle, Observations

The idea of a planetary origin for the solar cycle dates back to the nineteenth century. Despite unsurmounted problems, it is still advocated by some. Stefani, Giesecke, and Weier (2019) thus recently proposed a scenario based on this idea. A key argument they put forward is evidence that the ∼ 11 years-solar cycle is clocked. Their demonstration rests upon the computation of a ratio proposed by Dicke (1978) applied to the solar cycle time series of Schove (1955). I show that their demonstration is invalid, because the assumptions used by Schove to build his time series force a clocked behaviour. I also show that instabilities in a magnetized fluid can produce fluctuation time series that are close to being clocked.

Introduction

The recent opening and publication by Courtillot, Le Mouël and Lopes (Malburet, 2019) of a pli cacheté (sealed letter), entrusted to the French Academy of Science by Jean [START_REF] Malburet | Sur la période des maxima d'activité solaire[END_REF] [START_REF] Malburet | Sur la période des maxima d'activité solaire[END_REF], highlights the early interest in the search for a link between solar cycles and tides. In fact, the work of Malburet was already known from a detailed and luminous report he wrote (in french) for the journal 'L'Astronomie' in 1925 [START_REF] Malburet | Sur la cause de la périodicité des tâches solaires[END_REF] 1 . In that report, he correctly estimates the tidal forcing of planets on the Sun as being proportional to m p /d 3 p where m p and d p are the mass of a planet, and its distance from the Sun, respectively. This scaling yields tidal forcings proportional to ∼ 4.0, ∼ 3.8, ∼ 1.8, and ∼ 1.7, for Jupiter, Venus, Earth and Mercury, respectively (the contributions of the other planets are at least 10 times smaller). With this in mind, Malburet shows a good correspondance between the dates of solar maxima and the dates of 'weak deviations from Jupiter-Venus-Earth syzygies2 '. Malburet's idea was later taken up and extended by [START_REF] Wood | Physical sciences: sunspots and planets[END_REF] and [START_REF] Okhlopkov | The gravitational influence of Venus, the Earth, and Jupiter on the 11-year cycle of solar activity[END_REF].

There are however two serious problems with this idea:

1. The amplitude of the tidal forcing on the Sun is extremely small (< 10 -8 kg m -3 ), yielding accelerations one thousand times smaller than observed in the convective zone of the Sun (De [START_REF] Jager | Do planetary motions drive solar variability?[END_REF]. The forcing is 100 000 times smaller that the tidal forcing of the four Gallilean satellites on Jupiter (which is similar to the tidal forcing of the Moon on Earth). 2. The ∼ 11.2 years-period inferred from the 'weak deviations from Jupiter-Venus-Earth syzygies' is an artificial construction which has no signature in the complete tidal signal, as demonstrated by [START_REF] Okal | On the planetary theory of sunspots[END_REF], and illustrated in Appendix A.

It is therefore surprising that this idea receives a renewed attention [START_REF] Scafetta | Does the Sun work as a nuclear fusion amplifier of planetary tidal forcing? A proposal for a physical mechanism based on the mass-luminosity relation[END_REF][START_REF] Okhlopkov | The gravitational influence of Venus, the Earth, and Jupiter on the 11-year cycle of solar activity[END_REF][START_REF] Baidolda | Search for planetary influences on solar activity[END_REF][START_REF] Courtillot | On the prediction of solar cycles[END_REF][START_REF] Charbonneau | External Forcing of the Solar Dynamo[END_REF]. Even more surprising is the enthusiasm shown by Stefani and colleagues who have published no less than 7 articles exploiting this idea (see [START_REF] Stefani | Shaken and Stirred: When Bond Meets Suess-de Vries and Gnevyshev-Ohl[END_REF] and references therein). It seems that the main reasons that give confidence to these authors is their demonstration that the solar cycle is clocked, and probably the belief that this property requires a clocked forcing that only planetary motions can provide. Their demonstration, exposed in [START_REF] Nataf | A model of a tidally synchronized solar dynamo[END_REF], rests upon the computation of 'Dicke's ratio' of a 1000-years long time series of solar minima, which favors a clocked origin over a random-walk type origin.

The main objective of the present article is to show that the demonstration of [START_REF] Nataf | A model of a tidally synchronized solar dynamo[END_REF] is invalid. I also show examples of fluid instabilities that naturally produce clocked-looking time series.

The demonstration of Stefani et al (2019)

Stefani, Giesecke, and Weier (2019) picked up an idea that [START_REF] Dicke | Is there a chronometer hidden deep in the Sun?[END_REF] proposed for testing whether there is a "chronometer hidden deep in the Sun". The aim of Dicke was to distinguish a clocked behaviour from an 'eruption hypothesis', in which solar cycles would appear with a random phase. While Dicke restricted his analysis to the post-1705 time series of 25 solar maxima, [START_REF] Nataf | A model of a tidally synchronized solar dynamo[END_REF] extend it to 92 solar cycles starting in AD 1000, in an attempt to obtain a better statistical significance.

Dicke's ratio

Dicke (1978) noticed that a succession of 3 very short cycles starting in 1755 was followed by a very long cycle, as if the Sun was trying to keep up with some internal clock. He proposed several statistical tools to assess the existence of such a clock.

Consider a time series t(i) of N events i. In a perfectly clocked time series, all event dates t(i) are evenly spaced. When gaussian distributed noise is added, each event date is displaced from the regular grid by some random time, yielding a corresponding distribution of cycle durations d(i) = t(i) -t(i -1). In contrast, when events occur with a random phase, cycle durations d(i) have a gaussian distribution, and event dates are obtained as t(i) = t(i -1) + d(i). Clocked and random-walk time series can yield the same mean cycle duration d and standard deviation σ, but their statistical properties are not all identical. [START_REF] Dicke | Is there a chronometer hidden deep in the Sun?[END_REF] introduced a ratio that measures this difference, which Stefani, Giesecke, and Weier (2019) refer to as 'Dicke's ratio'. Dicke's ratio Di(n) computed for subsets of n ≤ N consecutive events is defined by:

Di(n) = n i=2 δ 2 n (i) n i=2 (δ n (i) -δ n (i -1)) 2 , (1)
where δ n (i) = t(i) -f n (i) is the deviation of event date t(i) from a best linear fit f n (i) = a n i + b n of the n dates.

According to [START_REF] Dicke | Is there a chronometer hidden deep in the Sun?[END_REF], the expectation of Dicke's ratio is:

E(Di clock (n)) = n 2 -1 2(n 2 + 2n + 3) -→ 1 2 when n → ∞ (2)
for a clocked time series, while it is:

E(Di rand (n)) = (n + 2)(n 2 -1) 3(5n 2 + 6n -3) -→ 1 15 n when n → ∞ (3)
for a random-walk time series. The expectation of Di(n) is independent of d and σ for both families, but the spread around the expectation does depend upon σ. [START_REF] Dicke | Is there a chronometer hidden deep in the Sun?[END_REF] applied this statistical tool to the time series of sunspot numbers starting in 1705. Due to the limited number of cycles, he could not reach a very definitive conclusion.

Schove's solar cycle time series

To reach a firmer conclusion, [START_REF] Nataf | A model of a tidally synchronized solar dynamo[END_REF] complemented the post-1705 solar minima series by solar minima dates as far back as AD 1000, following [START_REF] Schove | The sunspot cycle, 649 BC to AD 2000[END_REF]. Indeed, Schove published in 1955 the outcome of a very ambitious venture: dating maxima and minima of the solar cycle from 653 BC to AD 2025! Pre-1705 observations of sunspots being very rare, he mostly relied on reports of the observation of polar aurorae. In order to make up for the limited amount of reliable data, [START_REF] Schove | The sunspot cycle, 649 BC to AD 2000[END_REF] explicitly mentions (p.131) that he made use of two assumptions to build his 26-century-long table. These assumptions are reproduced in Figure 1. [START_REF] Schove | The sunspot cycle, 649 BC to AD 2000[END_REF] to construct his series of solar maxima, as displayed in his article p.131.

Arguments for a rebuttal

Schove's assumption (b) as listed in Figure 1 clearly suggests that his time series of solar maxima is clocked by construction. In order to be more specific, I have built synthetic solar cycles to test the impact of Schove's assumptions on the character of the resulting time series.

Synthetic solar cycles

The well-documented 24 solar cycles from 1755 yield a cycle duration (time between maxima) of 11.0 ± 2.0 years. Extending backwards to AD 1000 with Schove's dates yields 92 maxima separated by 11.1 ± 2.2 years. I have built two different families of synthetic cycles: a random-walk family, and a clocked family. Both retain the post-1755 dates of solar maxima, as distributed by WDC-SILSO, Royal Observatory of Belgium, Brussels (SILSO World Data Center, 1755-2019).

• The random-walk series are built by drawing at random normally distributed cycle durations, with a mean of 11.1 years and a standard deviation of 2.0 years. The dates of maxima are then constructed by cumulative difference from the date of the oldest post-1755 maximum. • The clocked series are built by extending the post-1755 dates of maxima backwards in time with a constant duration of 11.1 years, and then adding to the obtained dates a normal random noise with zero-mean and a standard deviation of 1.5 years, this value providing the desired standard deviation of 2.0 years for cycle durations.

Figure 2 displays the probability density function (pdf) obtained with 10 000 realizations, for both the random-walk series (Figure 2a) and the clocked series (Figure 2b). The pdf of Schove series is also drawn in Figure 2b.

The impact of Schove's assumptions

I have examined which of the 10 000 realizations of both families comply with the two assumptions of [START_REF] Schove | The sunspot cycle, 649 BC to AD 2000[END_REF] recalled in Figure 1. I find only 3 random-walk realizations, and 42 clocked realizations. In other words, the assumptions used by [START_REF] Schove | The sunspot cycle, 649 BC to AD 2000[END_REF] practically exclude random variations of the duration between solar maxima.

The cycle duration pdf of Schove-compliant realizations are shown in Figure 2, while their time series and deviations from a linear fit are displayed in Appendix B together with those of Schove's series. Figure 3 shows Dicke's ratios3 for a random selection of 100 realizations of both families. The thick black line shows the mean of Dicke's ratio for all realizations, while the thick red line is the mean of the Schove-compliant realizations. Note that Dicke's ratio of the rare Schove-compliant random-walk series is much lower that the mean of all realizations, and that individual realizations can plot far from the mean.

Quasi-clocked magnetohydrodynamic instabilities

The original goal of [START_REF] Dicke | Is there a chronometer hidden deep in the Sun?[END_REF] was to test the compatibility of solar cycle time series with an 'eruption hypothesis' expressed by [START_REF] Kiepenheuer | The Sun[END_REF] as 'each cycle represents an independent eruption of the Sun which takes about 11 years to die down'. Fluid dynamics and magnetohydrodynamics instabilities do not necessarily behave that way, even when strong turbulence is present. For example, quasi-periodic magnetic oscillations have been reported in the VKS dynamo experiment [START_REF] Berhanu | Bistability between a stationary and an oscillatory dynamo in a turbulent flow of liquid sodium[END_REF] at Reynolds numbers above 10 7 . Nonetheless, we have seen that Dicke's ratio yields a more stringent measure of the clocked behaviour of a time series than provided by visual inspection or pdf.

The Grenoble DTSΩ liquid sodium experiment exhibits magnetic fluctuations that are often quite regular [START_REF] Schmitt | Rotating spherical Couette flow in a dipolar magnetic field : Experimental study of magneto-inertial waves[END_REF]. From one such experiment, I have extracted time series of maximum induced magnetic intensity (see Appendix C for details), and computed Dicke's ratio of several sequences of 100 consecutive maxima. Figure 4 shows that the behaviour of these magnetic fluctuations is quasi-clocked, even though the Reynolds number is 8 × 10 6 and the standard deviation of the fluctuations about 30%.

Conclusion

The demonstration by [START_REF] Nataf | A model of a tidally synchronized solar dynamo[END_REF] of a clocked behaviour for solar cycles is invalid because the 1000 years long time series they use [START_REF] Schove | The sunspot cycle, 649 BC to AD 2000[END_REF] is clocked by construction. The astrological quest for a link between solar cycles and planetary tides remains as unfounded as ever. Magnetohydrodynamics instabilities can produce quasi-periodic fluctuations that appear as almost clocked. Tidal minima are better seen in Figure 5b, which displays the envelope of the tidal signal. The minimum of the series (at t = 1245 days) is marked by a blue dot, while the second maximum (at t = 8174 days) is shown by a red dot. The corresponding positions of the 4 planets are given in Figures 5c and5d; respectively. As expected, all 4 planets are almost aligned with the Sun at the maximum, while Jupiter cancels Venus' tide and Mercury cancels the Earth's tide at the minimum.

Figure 6 presents the spectra of a 360 years-long synthetic tidal signal (blue) and of its envelope (orange). The 4 peaks of the former simply correspond to half a solar day as seen from each planet (see equation 4). The spectrum of the envelope is dominated by peaks at the periods of syzygies of pairs of planets with the Sun, and their overtones. The spectrum is almost flat for periods beyond 300 days. This plot mimics Figure 3 of [START_REF] Okal | On the planetary theory of sunspots[END_REF], which shows the full tidal spectrum, 'taking into account the complete orbital elements [of the 4 planets], including eccentricity, inclination and their variation with time' over a period of 1 800 years. Orbital eccentricity adds up small tidal peaks at the orbital periods of Jupiter (11.86 years) and Mercury (0.241 years), but nothing shows up at 11.2 years (4088 days), as emphasized by [START_REF] Okal | On the planetary theory of sunspots[END_REF]. . Power spectra of the tidal signal (blue) and of its envelope (orange), as a function of period. Tidal peaks of the envelope spectrum appear at syzygies of pairs of planets with the Sun. Following [START_REF] Okal | On the planetary theory of sunspots[END_REF], I label them with the initials of the two planets. The largest tide is 'VJ' at 118.5 days when Venus and Jupiter are aligned with the Sun.

B. Time series and deviations of synthetic solar maxima

Figures 7 and 8 display the time series and deviations of the synthetic solar cycles that comply with Schove's assumptions, to compare with Figure 1 of [START_REF] Nataf | A model of a tidally synchronized solar dynamo[END_REF]. Deviations of each realization are the difference between dates of maxima and a linear fit of these dates.

C. Dicke's ratio of magnetic fluctuations in the DTSΩ experiment

The DTSΩ experiment was built to study magnetohydrodynamics in the magnetostrophic regime, in which Lorentz and Coriolis forces are dominant. Fifty liters of liquid sodium are enclosed in a spherical container that can rotate around a vertical axis. An inner central sphere can rotate independently around the same axis, and hosts a strong permanent magnet producing an axial dipolar magnetic field. The 3 components of the induced magnetic field are measured at the surface of the outer sphere at 20 equally-spaced latitudes from -57 • to 57 • (see [START_REF] Schmitt | Magneto-Coriolis waves in a spherical Couette flow experiment[END_REF] for more details). The frequency spectra of electric and magnetic fluctuations reveal a quasi-periodic behaviour that can be linked to the presence of magneto-Coriolis waves [START_REF] Schmitt | Rotating spherical Couette flow in a dipolar magnetic field : Experimental study of magneto-inertial waves[END_REF][START_REF] Schmitt | Magneto-Coriolis waves in a spherical Couette flow experiment[END_REF] or instabilities [START_REF] Figueroa | Modes and instabilities in magnetized spherical Couette flow[END_REF][START_REF] Kaplan | Dynamic domains of the Derviche Tourneur sodium experiment: Simulations of a spherical magnetized Couette flow[END_REF].

An example of such quasi-periodic magnetic fluctuations is given in Figure 9. The spin rate of the outer sphere is f o 10 Hz, while the differential rotation of the inner sphere is ∆f -20 Hz, yielding fluid velocities U 25 m/s. With an outer radius r o = 0.21 m, the Reynolds number Re = U ro ν 8 × 10 6 and the magnetic Reynolds number Rm = U ro η 60, with ν the kinematic viscosity, and η the magnetic diffusivity. Figure 9a displays a latitude-versus-time colorcoded image of the squared azimuthal magnetic fluctuations averaged over one turn of the inner sphere, in a time-window of some 250 rotations of the inner sphere. Latitudinally-coherent quasi-periodic fluctuations of variable intensity dominate. Figure 9b is the time-record of the same fluctuations at a latitude of 9 • . I have extracted the 100 maxima of this record, plotted as triangles in figure 9b. The duration between maxima is 0.11 ± 0.03 s. Dicke's ratio of the resulting time series is plotted in Figure 4, together with Dicke's ratio of the previous and next series of 100 maxima. All three series appear to be closer to a 'clocked' behaviour than to a random-walk behaviour.

Figure 1 .

 1 Figure1. The two assumptions made by[START_REF] Schove | The sunspot cycle, 649 BC to AD 2000[END_REF] to construct his series of solar maxima, as displayed in his article p.131.

Figure 2 .

 2 Figure 2. Probability density function of the duration of synthetic solar cycles. (a) Random-walk synthetics: gaussian distribution of all 10 000 realizations (teal), and pdf of the 3 realizations that comply with Schove's assumptions (magenta). (b) Clocked synthetics: nearly gaussian distribution of the 10 000 realizations (teal), pdf of the 42 Schove-compliant realizations (magenta), and pdf of Schove cycle durations (red).

Figure 3 .

 3 Figure 3. Dicke's ratios of (a) random-walk and (b) clocked synthetic solar cycles. Random selection of 100 realizations (grey), mean of all 10 000 realizations (thick black), Schove-compliant realizations (red) and their mean (thick red), Schove time series (blue squares). The thick green and magenta lines show the expectation of Dicke's ratio for a random-walk (equation 3) and for a clocked law (equation 2), respectively. Note that all synthetics share the same post-1755 reliable time series.

Figure 5 .

 5 Figure 5. Time series and envelope of tidal forcing as a function of time in a realistic synthetic solar system. (a) Tidal forcing at a given meridian on the Sun, calculated from equation 4 over 10 000 days (27.4 years). The periods of the 4 planets are marked by vertical lines of different colors. The maximum tide (achieved at t = 0) is shown by a horizontal dashed black line. (b) envelope of the tidal forcing. The minimum and (second) maximum are marked by a red dot and a blue dot, respectively. The periods of Jupiter and Venus are shown by vertical lines. (c) Positions of the 4 planets at the minimum tide (t = 1245 days). (d) Positions of the 4 planets at the (second) maximum tide (t = 8174 days).

  Figure6. Power spectra of the tidal signal (blue) and of its envelope (orange), as a function of period. Tidal peaks of the envelope spectrum appear at syzygies of pairs of planets with the Sun. Following[START_REF] Okal | On the planetary theory of sunspots[END_REF], I label them with the initials of the two planets. The largest tide is 'VJ' at 118.5 days when Venus and Jupiter are aligned with the Sun.

Figure 7 .

 7 Figure 7. (a) Time series and (b) deviations of Schove-compliant random-walk synthetic solar cycles. The blue squares are for Schove's maxima before 1755, and from SILSO World Data Center (1755-2019) after 1755.

Figure 8 .

 8 Figure 8. (a) Time series and (b) deviations of Schove-compliant clocked synthetic solar cycles. The blue squares are for Schove's maxima before 1755, and from SILSO World Data Center (1755-2019) after 1755.

Figure 9 .

 9 Figure 9. Magnetic fluctuations in the DTSΩ experiment. (a) Color-coded image of the square of azimuthal magnetic field fluctuations, normalized by the square of the imposed magnetic field (%) at the surface of the outer sphere. The magnetic field is measured at the 20 latitudes which form the y-axis. The x-axis is time given in number of turns of the inner sphere with respect to the outer sphere. The records are averaged over one such differential turn. (b) Extraction of the same signals at a latitude of 9 • , showing a succession of peaks labelled by magenta triangles.

alignment of all these planets with the Sun.

SOLA: Solar_cycles.tex;[START_REF] Charbonneau | External Forcing of the Solar Dynamo[END_REF] 17:40; p. 3

SOLA: Solar_cycles.tex;[START_REF] Charbonneau | External Forcing of the Solar Dynamo[END_REF] 17:40; p. 4

Note that all Dicke's ratios are calculated backwards in time as in[START_REF] Nataf | A model of a tidally synchronized solar dynamo[END_REF], since more recent dates are considered more reliable.

SOLA: Solar_cycles.tex; 28 June 2022; 17:40; p. 7 Nataf H.-C.

Acknowledgments I thank André Giesecke and Frank Stefani for providing clarifications on their computation of Dicke's ratio. I thank my colleagues of the geodynamo team of ISTerre for useful discussions and encouragements, and an anonymous reviewer for helpful suggestions. This article is dedicated to Emile Okal and to the memory of Don L. Anderson. : Solar_cycles.tex;28 June 2022; 17:40; p. 6

SOLA

Materials Availability All matlab scripts and data used to produce the figures of this article are available as supplementary material.

Conflict of interest

The author declares that he has no conflicts of interest.

Appendix

A. Frequency spectra of synthetic tidal forcings I have built synthetic tidal forcings to illustrate the lack of evidence for a 11.2 years periodicity, as demonstrated by [START_REF] Okal | On the planetary theory of sunspots[END_REF]. I model the tidal forcing exerted by Jupiter, Venus, Earth and Mercury, assuming circular orbits, and starting at a time when all planets are aligned. Mass, distances and orbital periods are taken from https://nssdc.gsfc.nasa.gov/planetary/factsheet/. Tidal forcing of this simplified planetary system at a given meridian on the Sun can be expressed as:

where m p and d p are the mass and distance from the Sun of planet p, respectively. T p is its orbital period, while T S is the duration of a solar day, taken as 27 days.

Figure 5a shows the predicted tidal signal over a period of 10 000 days (27.4 years). The orbital period of Jupiter (11.86 years) is indicated by vertical red bars. The maximum tidal forcing achieved at t = 0 is marked by a horizontal dashed line. It can be seen that forcings almost as large occur many times during one orbital period of Jupiter, typically each time Jupiter and Venus are aligned with the Sun.