Structural impact of the start-up sequence on Pelton turbines lifetime: Analytical prediction and polynomial optimization

Hydropower plants play a crucial role in the power system facing ambitious renewable energy targets. Due to their inherent controllability, they are well suited to provide flexibility to the grid. However, an increased flexibility provision leads to a prolonged operation in off-design conditions and more frequent hydraulic transient, notoriously detrimental in terms of dynamic loads induced on the hydraulic machines. The aim of this paper is to provide a methodology to estimate the damage on Pelton turbines due to the accumulated mechanical stress during their start-up sequence. Transient sequences, as the back-to-back start-up of centrifugal pumps in a pumped-storage power plant, requires the hydraulic machine to operate in conditions featuring high mechanical torque values, leading to an increased structural damage with respect to the nominal operating point. In view of the above, the paper proposes a polynomial optimization algorithm to decrease the fatigue-related operational cost induced on the runner during a start-up sequence. The performance of the algorithm is assessed on a real test case. The presented algorithm identifies the admissible injector opening rate as the crucial parameter in the start-up computation, and shows the dependence of the optimal opening law on the admissible rate with the bucket resonance characteristics.


Research context
In the last decades, a substantial rearrangement of the electricity industry in a market regulated structure has taken place in response to the transition from a vertical integration supply to the current unbundled system.Within this context, system operators have to simultaneously deal with the electricity negotiations occurring in dedicated market places [1] at different time scales, but also with the intermittent nature of the increasingly diffused non-dispatchable renewable energy sources (RESs) to reach the energy transition targets.As a consequence, the energy equilibrium of the power system must be ensured by the re-allocation of parts of the power generation units included in the day-ahead schedule agreed by the electricity market.The establishment of the constrained generation schedule requires therefore crucial redispatching measures, and hydropower plants (HPPs) are widely recognized to be key assets to provide the required generation flexibility to power systems.The existing literature [2][3][4][5][6] extensively discussed the sharp increase of flexibility demand of HPPs and their consequent intensified exploitation in transient regimes.However, this * Corresponding author.E-mail addresses: aldo.alerci@epfl.ch(A.L.Alerci), elena.vagnoni@epfl.ch(E.Vagnoni), mario.paolone@epfl.ch(M.Paolone).inherent flexibility forces HPPs to undergo increasingly frequent transient operating conditions, such as start-up or shutdown sequences, that are acknowledged to significantly contribute to the fatigue induced structural damage experienced by the hydraulic machines in operation.The structural damage caused by transient operating conditions on hydraulic machines has been studied for different turbine types, but an accurate damage quantification is hard to be performed.Such estimations are generally accomplished by performing experimental tests on full-scale machines.However, such tests require an interruption of the HPP operation to perform dedicated experimental campaign.An alternative approach relies on the transposition of results obtained from reduced scale model tests [7].As a consequence, the quantification of the fatigue accumulated during a transient operation is today not fully assessed and the modeling and the optimization of such transient trajectories are therefore rarely available for hydraulic machines operations.

Literature survey
Since the increasing need of HPPs flexibility provision is a relatively recent topic, the amount of studies about the transient induced https://doi.org/10.1016/j.renene.2023

Number of buckets (-)
Pelton turbines, their dynamic behavior during transient regimes has been studied in [13], which analyzed the vibration measured by accelerometers situated in different locations of the studied two-runners Pelton unit bearings during start-up sequences.The analyzed signals showed higher vibration levels during the start-up phase compared to the ''speed-no-load'' regime as well as for different steady state load cases.Despite the fact that these vibration levels have not been detected on-board the runner, the preliminary analysis of the vibrations transmissibility from the buckets to the bearings ensures that the considered scenarios distinctly depict the moment characterized by the most critical dynamic loads.Andolfatto et al. [14] employed a coupling of computational fluid dynamics (CFD) simulations, used to compute the jet pressure profiles inside the Pelton turbine bucket, with a structural finite element method (FEM) analysis to estimate the stress mean value and the peak-to-peak amplitude characterizing the solicitation at the bucket root during a start-up sequence.A Pareto front approach was subsequently applied to compute an optimal time of startup sequence, highlighting a trade-off between a low starting time and a low Usage Factor.
In [15], a detailed analysis of the modal characteristics of four Pelton turbine prototypes is presented, with a focus on the identification of dominant modes of vibration of the runners dynamic response.The comparison between the performed experimental and numerical study indicates that the bucket-dominated modes are predominant in terms of vibration amplitude compared to the disk-dominated ones, especially for circumferential vibrations.In relation to the vibration level, a failure analysis of a Pelton turbine bucket [16] provides an insight on the different phenomena which are responsible for the detected vibrations during the unit operation, essentially induced by the bucket passing excitation, centrifugal forces and the dynamic amplification due to the runner modal behavior.As a result, the bucket may be considered the most solicited part of the runner.This hypothesis is confirmed in [14,17], where the application of strain gauges at the bucket root, i.e. the hotspot concerning mechanical stresses in Pelton turbines, allowed measuring the impact of different operating conditions on the turbine structural integrity.

Paper's contribution
Previous research activities on the transient operation of Pelton turbines and correlated damage have been mostly restricted to experimental or computationally expensive (CFD and FEM) investigations to relate Pelton turbines operations with the quantification of dynamic loads.Furthermore, the limited amount of studies on hydraulic turbine start-up optimization, along with the increasing interest due to the evolving power system requirements, gives relevance to the topic.In view of the current state of the art, the main objective of this paper is twofold, as listed below: • analytical prediction of the mechanical stresses induced at the bucket root of a Pelton turbine during transient operating conditions; • formulation of an optimization problem aiming at computing an optimized start-up sequence minimizing the strain energy dissipated during the whole sequence.
As illustrated in [18][19][20][21][22][23], a dissipated strain energy quantification is particularly suitable to evaluate the induced damage on the structure in terms of lifetime reduction.By leveraging this finding, the present research elucidates the rootcause of structural damages in Pelton turbine during transient operations, and correlates the critical operating conditions by providing a new approach to predict and minimize the fatigue-induced damage during transients such as in a start-up sequence.In this regard, an analytical modeling of the dynamic loads is implemented and integrated in an optimization algorithm involving fast-growing convexification methods which represents a leap for the monitoring and control strategies of Pelton turbines operation in HPPs.
For this purposes, Section 2 describes the problem statements and introduces the developed methodology.Section 3 presents the main steps of the methodology, namely: the structural and hydraulic simplified models and the optimization algorithm to extract the optimal start-up sequence.In Section 4, the considered case study and the outcomes of the work are shown, whereas a discussion about the developed methodology and the obtained results are given in Section 5.

Problem statement
To predict the structural damage experienced by Pelton turbines, this study proposes a framework to model the turbine structure's critical location in terms of dynamic loads as an analytically tractable problem.More specifically, the considered critical location is the root of the bucket that links it to the turbine disk.To obtain a damage chart, the strain energy [24] dissipated at the bucket root is mapped on a mesh of steady-state operating conditions, defined by the rotational speed and the injector opening ratio, covering the entire operational domain of the turbine.The analysis of the transient-related damage associated to the start-up sequence1 is carried out by employing a quasi-steadystate approach.The transient sequence is discretized in time steps, and at each time step the turbine is modeled as operating at a steadystate condition.Therefore, the whole strain energy dissipated during a transient sequence is calculated as the sum of the instantaneous damage contributions related to each steady-state operating condition encountered throughout the sequence.The damage prediction is used by an optimization algorithm to compute a start-up sequence adjusting the turbine control variables to minimize the structural damage cumulated between given initial and final operating conditions.Since the optimization algorithm requires the analytical description of the objective function, which is extracted from the dissipated strain energy chart, and given the complexity of the latter, a surface fitting approach based on the Multivariate Adaptive Regression Splines (MARS) [25] is adopted.The treatment of the non-convex nature of this polynomial optimization problem (POP), notoriously NP-hard in general [26], necessitates a specific procedure.In this respect, a semidefinite relaxation is applied to tackle the non-convexity of the problem and, therefore, to identify the global optimal solution.

The bucket
The modeling of the mechanical stresses experienced by a solid structure under a specific load case requires a reduced-order approach to represent geometrical and structural properties.In the case of the Pelton turbine, the buckets configuration may be approximated as a series of cantilever beams, supported by the turbine disk at one end and free at the opposite.This simplified representation is widely employed in structural mechanics research concerning single-end-supported solid bodies.Typical examples are investigations on the structural behavior of aircraft wings, frequently reported in research papers and technical reports [27][28][29][30].
Given the preminence of bucket-dominated modes compared to the disk-dominated ones in terms of vibration amplitude [15], the disk-bucket system shown in Fig. 1 is modeled as a cantilever beam supported by an infinitely rigid wall as in Fig. 2. The rectangular-shaped constant section of the beam,   , is defined by the thickness of the turbine disk,   and by the circumferential segment corresponding to one bucket   =  ⋅   ∕  .Due to the non-uniform distribution of the bucket mass  in radial direction, the bucket is represented by a concentrated mass at the radial coordinate  = ∕2 and by a beam of negligible mass   ≪ , similarly to what is proposed by Schmied et al. [31].
With this approach, it is possible to resolve the bucket dynamic problem by considering the structure as a damped oscillator characterized by the following linear differential equation: where () is the mass displacement along the considered coordinate and  () is the external excitation term.Such computation implies two important limitations with respect to the physics of the problem: (i) with the given simplified representation of the bucket as a bending beam, (ii) by considering only the first bending mode, we reduce the complexity of the vibration assessment compared to the actual response.Nevertheless, given the shape and the orientation of the water jet excitation, the dynamic response of the bucket can be reasonably reduced to this mode [32].

Physical excitation
The excitation modeling is developed by considering the interaction between one bucket and a single injector.The excitation acting on the runner during an entire rotation is then computed including the amount and the location of the injectors of the considered Pelton turbine.As shown in Figs. 1 and 2, the excitation on the bucket can be divided in two components: the centrifugal force and the impinging jet force.In the simplified model, both are considered to act on the concentrated mass.In fact, although the location where the water jet impacts on the bucket depends on the bucket angular position (), the variation of the impinging angle between the jet and the bucket root surface causes a bending torque that can be approximated as a tangential force constantly applied at the radial coordinate  = ∕2.The detailed procedure to compute the jet force   () is provided in Appendix A.

Fourier decomposition of the excitation
To compute the dynamic response of the bucket in a steady-state operating condition, the different components of the external excitation are treated separately.The response related to the centrifugal effect of the excitation simply relies on the static structural behavior of the bucket.Instead, the jet contribution of the excitation must be characterized dynamically.Moreover, to compute the response of the bucket in bending mode, the dynamic amplification given by the interaction between the jet force and the modal characteristics of the bucket must be considered.Therefore, the jet force is decomposed in Fourier series to enable the detection of the potential match between the natural frequencies of the bucket and the different sinusoidal components of the excitation.Nevertheless, the periodicity of the Fourier series decomposition intrinsically weakens the influence of transient effects related to the structural damping.The computed dynamic response, especially during the free oscillation of the bucket between two subsequent jet impacts, may therefore be less reliable for lower values of rotational speed.
The Fourier analysis states that any periodic signal  () of period  = 2   can be decomposed in Fourier series terms [33] as follows: with 1 2  0 the static component corresponding to the mean value of the signal.The Fourier coefficients   and   are extracted respectively from the real and the imaginary part of the Fourier transform of the signal.The discrete Fourier series of the steady-state periodic force   () is numerically computed with a finite number  of harmonics.Since the position of the injectors may not be symmetric with respect to the runner axis, the period  of the excitation is defined by an entire rotation of the runner.

Steady-state dynamic response in bending mode
The steady-state solution of the bucket displacement in bending mode is computed by solving Eqs. ( 1) and ( 2).The linearity of Eq. ( 1) allows summing the solutions corresponding to each term of Eq. ( 2), leading to the following expression [33] of the displacement response: where the sinus and cosinus terms are grouped to introduce the −th harmonic amplitude   and the phase shifts   and   : (4)

Strain power density chart
The damage associated to the steady-state dynamic response is evaluated by computing the strain energy dissipation due to the stress at the root of the bucket.At this location, the reference section is given by   and   .The details of the calculation of the equivalent stress   are provided in Appendix B.
Once   is estimated, the associated strain energy density can be computed at any location of   .In this paper, the damage prediction is performed by considering the strain energy density at   =   2 and approximating the shear modulus value as  ≈  3 .The strain energy density can be computed [34] as: For a steady-state operating conditions, the equivalent tensile stress signal   () computed over a time lapse  − is used to derive the amount of strain energy density transferred to the considered location during this time lapse.The time average of the latter leads to the definition of the strain power density  − [W/m 3 ], which is defined as: The equivalent discretized formulation is the following: Since the period of the excitation  depends on the rotational velocity of the runner, the time lapse  − is adjusted to have the same integer number of runner rotations for every operating condition.
Given that the response sampling frequency is kept constant for each operating condition, the amount of samples  therefore varies between conditions featuring different values of .
To represent the structural damage as a function of the operating condition, the operating range defined by the runner rotational speed  and the injector opening   is discretized between the turbine standstill [ = 0 ;   = 0] and a regime exceeding the nominal rotational speed of 25% and featuring the complete opening of the injector,   = 1.The strain power density computed for these discrete time steps generates the damage chart shown in Fig. 3. Two main aspects can be commented: first, the general trend of the mechanical stress induced by the jet force and the centrifugal effect is clearly recognizable.If considered separately (i.e. with one of the two effect set at zero), each force induces a monotonous increase of the stress and, therefore, of  − .By contrast, the combination of non-zero  and   values produces a curved development of the stress, suggesting the existence of an optimal trajectory of the turbine start-up sequence that may minimize the cumulated stress between the standstill and the nominal condition (yellow mark in Fig. 3).Second, the presence of local spikes in the  − magnitude is also remarkable.These spikes are related to the resonance conditions of the bucket bending modes and they are linked to specific values of  which determines the fundamental frequency of the bending excitation.Moreover, the spikes are characterized by an oscillating behavior along   .This trend is the consequence of the variation in the excitation signal shape depending on the   value: in fact, the presence of specific terms in the Fourier decomposition that match with the bucket resonance frequency results from the shape of the decomposed signal.
Once the damage chart is defined, the structural damage induced on the bucket during a transient sequence may be estimated.This estimation is performed by considering the transient-induced damage as a sum of the infinitesimal contributions given by the steady-state operating conditions attained throughout the whole sequence.In the case of a start-up sequence, the amount of strain energy cumulated during the sequence time duration   may be expressed as follows: This quantity can be approximated as: with   =    .

Non-convex polynomial optimization problem 3.4.1. Multivariate regression model
The developed methodology seeks a quantification of the turbine damage in steady-state and transient operating conditions.The accomplishment of this objective can be exploited by an optimization problem aiming at computing a start-up sequence minimizing the dynamic loads experienced by the runner.For this purpose, an analytical expression of the strain power density chart and the optimization methods must be determined.Given the complexity and the non-convex nature of the strain power density chart on the operating range of the Pelton turbine, an approach based on the Multivariate Adaptive Regression Splines is selected to determine the analytical expression of the function  − (,   ).The amount of generated hinge functions combined with the highest polynomial degree of the formulation defines the accuracy of the fit.

Objective function
The objective function to be minimized must express the structural damage, in terms of strain energy density, cumulated during the sequence.Using Eq. ( 9), the objective function is formulated as follows: with (  (  ) ,   (  ) ) the operating condition of the turbine at the −th time step represented by the decision variables of the problem,  and   .

Mechanical system constraints
To represent the turbine start-up sequence, the problem must be properly constrained.This implies the definition of initial and final operating conditions, as well as the turbine motion law and the physical limits of the system.Each constraint must be expressed as a function of the decision variables.The time duration of the start-up sequence   is set as an imposed parameter of the problem.
The initial and final conditions imposed on  are defined by: ≈ 0 and (   ) ≥   (11) whereas regarding   , the following conditions apply: The exact final state of the turbine rotational speed (   ), being determined by the motion law, is not set to a numerical value to bypass the risk of an unfeasible solution.A nearly-constant final rotational speed and a lower bound   are set.Since this paper investigates the acceleration phase of the turbine, without considering the synchronization to the grid, it has been decided to constrain the injector final state to its nominal value,  , , rather than the final rotational speed value set by the grid frequency.However, with an accurate representation of the motion law, (   ) will obviously be very close to the nominal rotational speed.The interdependence of the decision variables is revealed by the law of rotating motion: with  the rotational inertia of the system composed by the Pelton turbine and the generator rotor [35] and the term   ⋅  representing a load transmitted by the turbine, e.g. in the back-to-back start-up procedure of a multistage pump in a pumped-storage HPP.The term   includes the total number of activated injectors.The differential Eq. ( 13) must be expressed as a function of the decision variables and discretized by a finite differential method in a series of algebraic constraints applied to each time step.By expressing   as a function of the discharge (see Appendix A), and performing a second degree polynomial fit (  ()) on the characteristic curve of the turbine, Eq. ( 13) may be entirely written as a function of  and   , as follows: The variation of the injector opening rate value   with time is determined by both the physical reactivity of the needle's movement and the safety limits imposed by the HPP operators to avoid detrimental consequences induced by hydraulic transient phenomena.This constraint is formulated as follows: with   a constant value.
Regarding the instantaneous mechanical stress produced on the bucket root, a constraint on the equivalent tensile stress represented by   based on the yield strength of the considered material may be suitable.However, during the earlier computation of the steady-state damages, the highest peak values detected for   attained about one order of magnitude lower than the yield strength of the considered steel.For this reason, this constraint has been neglected.

Polynomial optimization problem
The polynomial optimization problem is defined by the minimization of the objective function (10) and by the constraints presented in Eq. ( 11)- (15), and it is expressed as follows:

Solution methods
The non-convex nature of the objective function, resulting from the interaction between the dynamic response of the bucket and its resonance conditions, prevents the computation of a global optimal solution using first order iterative algorithms as the gradient descent method.As a result of the  − fit performed by the MARS technique, the optimization problem ( 16) may be classified as a non-convex POP, NP-hard to solve in general.To treat such non-linear optimization problems described by polynomials several non-linear programming (NLP) methods exist [36], considering also integer variables in the case of mixed-integer non-linear programming (MINLP) methods [37] to be implemented for use in a branch-and-bound framework, or based on generalized reduced gradient (GRG) methods [38].However, checking the convergence of such algorithms to the global optimum for non-convex problems remains a difficult task, and the quality of the obtained solution can hardly be evaluated.By contrast, fast-growing methods relying on semidefinite relaxations, i.e. introducing positive semidefinite matrices in the constraints of the problem, allow for the computation of the global optimal solution under specific conditions.In this regard, Lasserre provides a methodology based on the definition of a hierarchy of semidefinite relaxations whose optimal values form a monotone sequence of lower bounds converging to the global optimum also in non-convex problems [39].Therefore, the use of this relaxation technique on the optimization problem presented in this paper guarantees that the solution provided by the optimization algorithm has converged to the global optimum.
The main idea is to reformulate the non-linear minimization problem as a linear program, hence convex, via semidefinite relaxations.The initial problem is converted in a primal semidefinite program (SDP) and in its dual problem.The primal problem is defined via the Riesz linear functional   ∶ R[] → R and it is constrained by moment and localizing matrices semidefinite positive.To implement the SDP relaxations, few algorithms exist in the literature.One of them, called GloptiPoly and developed by Henrion and Lasserre [40], relies on this approach.

MARS model constraints
The application of the SDP relaxation method to the initial minimization problem requires the reformulation of the latter in matrix form.However, the objective function   does not have the form of a standard polynomial expression: in fact, the MARS formulation of  − introduces the hinge functions form  (0,   − ) in the polynomial terms, with   the physical decision variable of interest and  the constant threshold of the considered hinge function.To explicit   in matrix form, each hinge function must be reformulated.The use of binary variables   is introduced to rewrite the hinge functions form as follows: for  = 1, … , nr HF (17) with h.v. an arbitrary high value, at least greater than the upper bound of (  − ), and nr HF the amount of hinge functions contained in the MARS formulation.

Relaxed SDP
Once the binary variables introduced by the MARS model are specified, the set of constraints of the POP is entirely defined.The number of constraints included in ( 16) can be expressed as follows: = 5 + 2 ⋅ (  − 1) + 3 ⋅ nr HF ⋅   (18) and the number of variables of the problem reads as: The relaxation procedure illustrated by Lasserre [39] is employed to relax (16) (21) for finitely many real non-zero coefficients   which are, as well as the degree   of the objective function, dependent on the selection of the MARS formulation employed to fit the  − function.
To write ( 16) as a relaxed SDP relaxed, let us introduce the vector   () defined as: Let us group the inequality constraints in the set of real-valued polynomials   () ≥ 0 and the equality constraints in the set of realvalued polynomials ℎ  () = 0, and introduce the sum of squares (SOS) decomposition to rewrite them in terms of semidefinite matrices.The SOS decomposition allows writing a generic polynomial () ∈ R[] 2 by means of   () and a positive semidefinite matrix  ∈ R ()×() such that: with () ∶= (   +
To formulate the relaxed SDP, the introduction of the Riesz linear functional   ∶ R[] → R is required to rewrite a non linear polynomial as a linear one.In the case of the objective function  (), each monomial   is replaced by a single real variable   as follows: Hence,   ( ) is expressed as a function of the real sequence  = (  ) ⊂ R and the same non-zero coefficients   that define the POP objective function.Furthermore,   allows defining the real symmetric moment matrix as: meaning that   is applied to each entry of the matrix   ()    ().
The equality and inequality constraints of (16), converted in to positive semidefinite matrices conditions in (24), are expressed in the relaxed SDP defining the localizing matrices, i.e. the moment matrices associated to the set of polynomial inequality constraints,  −  (  ), respectively to the set of equality constraints,  −  (ℎ  ).
The POP (16) can therefore be written in its relaxed form.The th Lasserre's SDP relaxation is expressed as: with  = max(  ∕2,   ,   ) the relaxation order of the polynomial SDP.
The sequence  is a feasible solution of the relaxed SDP problem (27), with the value of the relaxed objective function   ( ) =  ().The inequality   ≤  * holds for every .Moreover, since an increasing  leads to a more constrained SDP, the relation  +1 ≥   holds as well.Hence, (27) defines a hierarchy of semidefinite relaxations of the initial problem (16).For non-convex problems, the finite convergence provided by the Putinar Positivstellensatz is generic, i.e. is ensured for a relaxation order  → ∞.Unfortunately, the size of the SDP relaxations increases rapidly with the size of the original problem, for an increasing relaxation order.This size explosion issue is due to the semidefinite matrices introduced in the relaxed problem.To handle problems characterized by a large size, the sparsity of matrices can be exploited.In this regard, Waki et al. propose the algorithm SparsePOP [41], a Matlab implementation of the sparse SDP relaxation method running by default the primal-dual interior-point solver SeDuMi (Self-Dual-Minimization) [42].Since the methodology employed in this paper requires an algorithm able to work with   and   values on the order of 10 3 , this problem size can be handled uniquely by exploiting sparsity.Therefore, the SparsePOP algorithm is selected.

HPP characteristics
The case study considered in this paper is one of the eight identical Pelton turbines installed in Veytaux I, powerhouse of the Forces Motrices Hongrin-Léman (FMHL) pumped-storage power plant located in the western region of Switzerland [43].This HPP features two powerhouses, for an overall installed power of 480 MW.In the selected powerhouse Veytaux I covering half of the total installed power, there are four horizontal axis ternary units, each one equipped with a centrifugal pump and two Pelton runners.Every Pelton runner features a rated power of 30 MW.In generating mode, an average gross head of 853 mWC for a nominal discharge of approximately 4 m 3 s −1 per runner are considered.The studied Pelton runner features two injectors at 80 degree from each other, and operates at a nominal rotational speed of 600 min −1 .The losses throughout the hydraulic circuit up to the injectors are represented by the coefficient  ℎ , estimated by considering a 2.5% loss at the rated power.
To define the structural characteristics of the runner, some approximations are required.The physical properties of the steel are considered.The available 2D and 3D geometries allow computing the mass of a single bucket and of the entire runner.The damping ratio  is estimated assuming the value of the relative damping  = 0.05% [31].
To model the rotational dynamic of the runner, the statistical analysis performed on the nominal transmitted torque as a function of the generator mass illustrated in [35] is adopted to estimate the latter, assumed as approximately 50 × 10 3 kg.The rotational friction coefficient   is estimated through the real-time data extracted from the Supervisory Control And Data Acquisition (SCADA) of the HPP.Finally, the coefficient   is determined by considering a load that increases linearly with the rotational speed until the rated power at the nominal operating condition is transmitted by the turbine.

MARS fit reduction
The global optimal start-up trajectory is provided by the SDP relaxation solution computed on a MARS-fitted surface which should be the closest possible to the strain power density chart shown in Fig. 3.This solution is ensured to converge to the global optimum for  → ∞.However, to have a problem computationally treatable a trade-off must be considered between (i) the refinement and accuracy of the MARS fit and (ii) the relaxation order  which must be limited.These parameters should be fine-tuned to obtain a satisfying quality of the SDP relaxed solution in a reasonable computational time.It is worth mentioning that the relaxation order  is strictly related to the degree   and   = 2 ⋅ [  ,   ] of the polynomials describing respectively the objective function  and the constraints of the initial problem by the following relation [39]: As a result, a decrease of the MARS fit refinement also requires a lower relaxation order .Moreover, the introduction of   leads to double the objective function degree   compared to the MARS formulation degree.This aspect, therefore, plays a role on the relaxation order selection.
The solution computed via the SDP relaxation represents the global optimum solution associated to a strain power density chart reproducing only partially the original chart.This outcome is employed to initialize a gradient descent algorithm which refines the solution.Thereby, the dominant effects of the non-convex nature of the problem can still be tackled, and the final solution being close to the global optimum can be considered as a best-effort optimum.

Relaxed problem optimum
The MARS formulation selected to solve the SDP relaxed problem features 31 hinge functions, combined together to generate a thirddegree objective function.The fit function features a RSQ value of 99.93% and a relative error of 16.7% measured at the maximal absolute error location.The corresponding strain power density chart is shown in Fig. 4. Since this accuracy is strongly affected by the non-convex regions of the surface, the rotational speed operating range is restricted to ignore the large peaks observed in Fig. 3 for  > 65 rad/s.However, the excluded operating range is never attained during the start-up procedure of the Pelton turbine.
The crucial prerequisite that motivates the choice of this formulation is represented by the presence of the most relevant resonanceinduced peak observed between  = 30 rad/s and  = 40 rad/s.For lower accuracy formulations, the shape and the magnitude of the peak are too far from the original surface illustrated in Fig. 3, by inducing an erroneous weight of the non-convexity of the function.However, with such fitting quality, the relative amplitude of the peak is independent on   , and its oscillating behavior is not reproduced: this feature observed in the original database is tackled in a second phase of the modeling through the gradient descent implementation.
For the sake of comprehensiveness, at first only one start-up trajectory is presented.The SparsePOP algorithm is initialized with the parameters illustrated in Table 1.The physical parameters are set in order to both fulfill the operating characteristics of the full-scale turbine and to reproduce a realistic start-up procedure of the HPP.The value of   is defined such that the minimal time duration required to ensure a complete opening or closing procedure of the injector in the range   = [0, 1] would be 33 s.The value imposed to   aims to include a load carried by the turbine during the acceleration phase and that stabilizes at 28.2 MW at the nominal operating condition.This scenario may correspond to a back-to-back start-up sequence of a centrifugal pump in a ternary group featuring a rated power of 30 MW.The a priori choice of the gross head as being the average head of the HPP implies that the delivered power at the imposed   value and the nominal rotational speed is slightly lower than the rated power.
The computed start-up sequence is presented in Fig. 5.The trend of the decision variables  and   is plotted as a function of time in Fig. 5(a), and the induced damage associated to each time step of the start-up sequence is illustrated in Fig. 5(b).
The result highlights that by providing an a priori value for   , the algorithm imposes an initial stall in the low-damaging operating range and then set an acceleration regime featuring the maximal admissible opening rate of the injector.The initial stall condition lasts up to the instant the algorithm must enforce the turbine acceleration to satisfy  1.
the constrained operating condition at the final instant    , i.e. about 13 s in the simulated scenario.
During the acceleration, the operating range, characterized by the resonance-induced peak, is crossed at the maximum injector opening rate.The operating condition at the final instant fulfills the imposed constraints: the injector opening is adjusted to ensure a nearly constant rotational speed, that varies less than 0.2% at the last time step    , and that is in agreement with the settled value of   .
The influence that the injector opening rate has on the optimal start-up sequence is investigated by comparing the turbine acceleration regimes under varying   values.In Figs. 6 and 7, three start-up sequences computed for injector opening rate in the range   = [2 × 10 −2 , 4×10 −2 ] s −1 are presented, plotted against time and on the damage chart.The same characteristics of the previously analyzed accelerating behavior are clearly observed for each opening case.The time duration of the initial stall is consistently related to the value of   .The higher the admissible opening rate is, the shorter the required time for the acceleration regime can be, decreasing the time spent by the runner on the most relevant operating ranges in terms of structural damage.The trend of the structural damage   against the injector opening rate   is shown in Fig. 8.It is interesting to notice how a shorter start-up leads to a decrease of the damage, despite the fact that in sharper start-ups the turbine explores regions of the  − surface characterized by higher instantaneous damage values.
In case higher   values are considered, it is worth computing the optimal start-up sequence for a shorter time duration of the transient sequence.The investigated start-up regimes covers the range defined by   = [0.8× 10 −1 , 1.6 × 10 −1 ] s −1 , i.e. increasing by a factor 4 the previously studied opening rates, with a fixed time duration of   = 15 s and a discretization set with  = 0.3 s.It is important to remark that such high   values describe a set of acceleration regimes going beyond the admissible scenarios for an HPP similar in size to the one considered as case study.The time duration required to fully open or close the injector associated to these   are comprised between 12.50 and 6.25 s.However, relevant properties of the optimization algorithm that is implemented in the framework of this investigation may be distinctly highlighted for such shortened time ranges, from the point of view of both the relaxed solution and the gradient descent refinement.Figs. 9 and 10 illustrate how the opening law is adapted to fulfill the final operating condition constraint in the imposed time duration.The initial stall situation is visible similarly to the low   regimes, and an overshoot of   above its nominal value during the last 3 s is needed to provide the required final speed value.The longer the initial stall lasts, the higher the overshoot amplitude must be.The influence of the injector opening rate on the induced structural damage is presented in Fig. 11.

Sensitivity analysis on the solution feasibility
The convergence of the SDP relaxed problem solution to a feasible optimum strongly depends on the selected parameters.The scenario described in Table 1 is selected as reference for this analysis.On one hand, the tuning of the physical parameters allows fulfilling the constraints that define the initial and final operating conditions of the start-up sequence.This is the straightforward result of the runner's law of motion described in Eq. (13).By contrast, the simulation parameters have a relevant influence on the feasibility of the solution with respect to the MARS model constraints, Eq. ( 17).The relaxation order  and the time discretization  have a substantial impact on the value of the binary variables introduced by the hinge functions.Since  must be bounded for computational cost reasons, the accuracy in the computation of the binary variables   is limited.This means that when one of the decision variables  or   approaches the threshold value  of a specific hinge function, the associated binary variable   deviates from its current value, reaching the opposite value solely once the threshold value  has been crossed.Therefore, this procedure lasts more than a single time step.Fig. 12 illustrates this trend: the step from 0 to 1 and viceversa is not instantaneous.For the vast majority of the time steps, this inaccuracy leads to a deviation in the corresponding value of the computed objective function from its ideal value calculated with the   variables perfectly fulfilling their binary constraint.It has been noticed that a decrease in the selected  causes an increase in the error computed on the quantity   .However, a too coarse time discretization would lead to an unreliable calculation of   as well: a particularly delicate circumstance could be the insensitivity to the resonance-induced peak, inducing a loss of information regarding the non-convexity of the  − surface.Therefore, a fine-tuning procedure must be performed to choose the appropriate .
The selection of the suitable relaxation order  is performed exploring the root mean square error (RMSE) calculated between the   variables fulfilling perfectly the binary constraint and the   values extracted by the SDP relaxed solution.Fig. 13 illustrates this analysis.The RMSE is constant for orders  ≤ 3. Crossed this threshold, increasing the relaxation order leads to a decrease in the RMSE.This trend is consistent with the property of the Lasserre SDP relaxation  order presented in Eq. ( 28).With a third degree MARS formulation, converted in a sixth degree function by the introduction of the binary variables, the threshold value for  to improve the accuracy of the relaxed solution is  = 3.The curve's slope decreases for higher , attaining RMSE = 7.5% for  = 6, and a decreasing error rate lower than 0.2% between  = 6 and  = 7.In view of a raise of approximately a factor 5 in the computational time between  = 6 and  = 7, the selected relaxation order is  = 6.

Gradient descent refinement
The optimal solution obtained from the SDP relaxed problem is used as input for a gradient descent refinement for two main reasons: (i) the relaxed problem is computed on a MARS fit with bounded nr HF and   values, i.e. with a limited fitting quality, and (ii) this solution involves the computation of binary variables   which induce an error on the calculation of the objective function   .Compared to the SparsePOP algorithm, the function employed to implement the gradient descent algorithm requires a significantly lower computational cost.Therefore, the limits imposed on the nr HF and   values are less restrictive, and the MARS formulation used to compute the objective function can reproduce more reliably the damage chart.Furthermore, the error obtained in the SDP solution mainly concerns the convex region of the  − surface at low  and   values.These considerations suggest that the use of a gradient descent method may improve the obtained solution, by approaching more closely the global optimal solution of the original problem.Since the gradient descent method does not guarantee the identification of the global optimum, the final solution can be considered as a best-effort optimal solution.
The gradient descent algorithm is initialized with the interior-point method to converge towards the solution.The objective function is computed on the  − surface shown in Fig. 14, defined with a MARS fit featuring the values nr HF = 68 and   = 7.The fit function features a RSQ value of 99.98% and a relative error of 12.9% measured at the maximal absolute error location.With this formulation the oscillating behavior of the resonance-induced peaks is recognizable.This advantage is of major interest for the start-ups with the higher admissible   values.In fact, a steeper    slope implies a greater flexibility of the injector opening trajectory, which can be adjusted to take advantage of the variable amplitude of the peaks to reduce the structural impact throughout the whole sequence.To improve the significance of the refinement, in particular regarding the sensitivity of the sequence to the resonance peaks, the time discretization is enhanced compared to     1. Fig. 13.RMSE computed on the whole set of binary variables associated to the start-up sequence of the reference case set with the parameters listed in Table 1, plotted against the relaxation order .The threshold at  = 3 is visible.the time step used to extract the SDP relaxed solution.The lower computational cost required by the refinement procedure allow adjusting the time step keeping the numerical effort reasonable.The employed time steps are halved:  = 0.5 s for the low   scenarios, respectively  = 0.15 s for the higher   values.For this purpose, the SDP solution is simply interpolated on the finer time steps mesh and implemented as input for the gradient descent approach.
The refinement is firstly applied to the set of sequences illustrated in Figs. 6 and 7.The sequence obtained from the SDP relaxed problem is used as initial point for the gradient descent.Figs. 15 and 16 illustrate how the algorithm improves the SDP relaxed solutions previously computed.The most remarkable improvements are related to the initial stall.This is consistent with the expectations, since the greatest error on   which affects the objective function in the SDP relaxed solution occurs at this operating range.The constraints applied to the mechanical system prevent any adjustment to the initial condition and the last 5 s of the sequence.With such admissible injector opening rates, the turbine does not have the required flexibility to bypass the high-stress operating conditions related to the resonance-induced peak.The influence of the gradient descent refinement on the variation of   with   is shown in Fig. 17   values, the deployment of the maximal opening rate throughout the entire acceleration phase is not necessarily the optimal operating procedure.Being highly responsive, the needle governor is able to adjust the injector opening sequence bypassing the extremely damaging operating range, still fulfilling the final condition constraints.This is notably visible in Fig. 19, where the first stage of the acceleration is the result of two main elements: i) avoid working at conditions characterized by a high   value and a low rotational speed, that would be encountered with a sharp initial acceleration phase, and (ii) to later bypass the impact of the resonance peak, exploiting the drop in the peak's amplitude in the range  = [35,39] rad/s and   = [0.45,0.60].
Fig. 20 and Table 3 show the variation of the damage value   induced by the refinement parameter.In reason of the moderate acceleration that characterizes the initial phase of the sequence for higher   values, the refined solution is rather close to the SDP relaxed solution.This leads to a lower damage variation provided by the gradient descent refinement solutions for increasing admissible injector opening rates.Improvements lower than 1% confirm that the SDP relaxed optimal solutions match rather accurately the damage values considered as the best-effort optimum computed by the refinement.
In view of the simulated regimes, the existence of a threshold value in the admissible opening rate   can be deduced: for   values higher than this threshold, the use of the maximal acceleration does not represent the optimal strategy.Considering the presented results, it is reasonable to consider that this threshold value is in the range   = [0.8× 10 −1 , 1.2 × 10 −1 ] s −1 .For start-up sequences characterized by   values higher than this threshold, the choice of   also has an impact.A smooth reduction of   from its optimal value implies an increase of the structural damage, until   reaches the lowest possible value which allow the turbine attaining the nominal operating condition.A further reduction would lead to an infeasible solution.This trend is illustrated in Table 4: the increase in the damage   with respect to the case   = 15 s is illustrated.The parameter   l.b. indicates the lower bound for the admissible opening rate to compute a feasible solution for the imposed   .

Conclusions
The paper presented a method to model and minimize the mechanical stresses induced at the bucket root of a Pelton turbine during transient operating conditions.The strain power density dissipated at this hotspot is computed on the turbine operating range by considering a mesh of steady-state operating conditions and decomposing the external excitation in Fourier series.Taking into account the limitations mentioned in Section 3.1 this approach allows for the extraction of the dynamic response in bending mode of the bucket and to map the structural damage on the entire operating range.
The damage chart is used to feed an optimization algorithm to minimize the cumulated structural damage on Pelton turbines during the start-up sequence.The non-convexity of the problem, introduced by the modal behavior of the turbine, requires the problem formulation to be relaxed to be tractable.The MARS regression analysis, applied to the strain power density mapping, and the Lasserre approach, based on a hierarchy of semidefinite relaxations, have been suitably integrated to formulate and solve the proposed non-convex problem.
The optimal start-up sequence is firstly computed with an objective function of bounded degree to be computationally tractable.This solution provides the estimation of the relaxed problem global optimum which is used as the initialization of the non-relaxed problem, solved via a gradient descent approach, to reliably reproduce the full modal properties of the turbine.Despite the fact that the gradient descent method does not ensure the convergence to the global optimum, this  limitation has been accepted to maintain a reasonable computational cost for the relaxed problem resolution.In case of availability of an important computing power, the full-degree objective function can be considered to calculate the global optimum of the relaxed problem, eliminating the need of the gradient descent refinement.
With the developed methodology, the influence of the HPP characteristics and of the runner structural properties can be identified and related to the control parameters employed to define the start-up sequence trajectory.In particular, it has been shown that the optimal injector opening procedure strongly relies on the admissible opening rate allowed on the HPP and on the resonance properties of the bucket.For low admissible opening rates, the optimal strategy aims at reducing the time duration of the sequence at most, by decreasing the turbine exposition to the transient-induced mechanical stresses.By contrast, with admissible injector opening rates higher than a threshold value, the deployment of the maximal acceleration cannot be considered as the optimal solution.This distinction is justified by the highlydamaging regimes that could be attained by a sharp acceleration with such opening rates and by the flexibility that the injector can exploit to prevent the runner to cross operating conditions causing important resonance-induced stresses.It is interesting to notice that this behavior, observed for opening rates higher than the threshold, shows an inverse correlation between the start-up time duration and the cumulated damage as retrieved by Andolfatto et al. [14], who applied a Pareto front optimization to resolve the a-priori formulated trade-off between a shorter start-up time   and a lower structural damage.In this regard, it is important to remark that the threshold value of the injector opening rate is strongly related to the considered rotational inertia of the system  , which is essentially due the generator size, as shown in Eq. ( 13).In fact, a higher rotational inertia lowers the threshold value of   because of a lower angular acceleration of the runner, leading to the appearance of a threshold for   values higher than 15 s.In this paper, the size of the generator has been estimated considering the rated power of the turbine, i.e. 30 MW.However, this equivalence may not be always justified.For instance, in the Veytaux I HPP considered as case study, every generating unit features two runners sharing a generator of a rated power of 60 MW.The load carried by a single runner may therefore vary in case of an unequal power balance between the two runners.

Fig. 3 .
Fig. 3. 3D map (a) and orthogonal projection (b) of the strain power density chart as a function of  and   computed with the analytical methodology.The yellow mark represents the nominal operating condition.

Fig. 4 .
Fig. 4. 3D map (a) and orthogonal projection (b) of the strain power density chart as a function of  and   -MARS with   = 31,   = 3.The yellow mark represents the nominal operating condition.

Fig. 5 .
Fig. 5. SDP relaxed solution.Start-up sequence trajectory of the reference case, as a function of time (a) and on the damage chart (b), computed with the parameters listed in Table1.

Fig. 6 .
Fig. 6.SDP relaxed solution.Time evolution of the decision variables   (a) and  (b) for three regimes featuring different admissible injector opening rates   and   = 45 s.

Fig. 7 .Fig. 8 .
Fig. 7. SDP relaxed solution. and   start-up trajectory on the strain power density chart for three admissible injector opening rates   .

Fig. 9 .
Fig. 9. SDP relaxed solution.Time evolution of the decision variables   (a) and  (b) for three regimes featuring different admissible injector opening rates   and   = 15 s.

Fig. 10 .Fig. 11 .
Fig. 10.SDP relaxed solution. and   start-up trajectory on the strain power density chart for different admissible injector opening rates   .

Fig. 12 .
Fig. 12.Time evolution of the 31 binary variables associated to the start-up sequence of the reference case computed with the parameters listed in Table1.

Fig. 14 .
Fig. 14. 3D map (a) and orthogonal projection (b) of the strain power density chart as a function of  and   -MARS with   = 68,   = 7.The yellow mark represents the nominal operating condition.

2
, and Table 2 summarizes the induced relative change in the damage value   for the illustrated start-up scenarios.Given the short Tab.Relative improvement of the gradient descent refinement on the cumulated damage for   = [2 × 10 −2 , 4 × 10 −2 ]s -1 .The variation induced by the gradient descent refinement is lower for the regime featuring a short-lasting initial stall, i.e. at the lowest   value.  [s −1 ]  ,  − ,   Relative improvement of the gradient descent refinement on the cumulated damage for   = [0.8× 10 −1 , 1.6 × 10 −1 ] s −1

Fig. 15 .
Fig. 15.SDP relaxed solution vs gradient descent (gd) refinement.Time evolution of the decision variables   (a) and  (b) for three regimes featuring different admissible injector opening rates   and   = 45 s.

Fig. 16 . 4
Fig. 16.SDP relaxed solution vs gradient descent (gd) refinement. and   start-up trajectory on the strain power density chart for different admissible injector opening rate   .

Fig. 18 .
Fig. 18.SDP relaxed solution vs gradient descent (gd) refinement.Time evolution of the decision variables   (a) and  (b) for three regimes featuring different admissible injector opening rates   and   = 15 s.

Fig. 19 .
Fig. 19.SDP relaxed solution vs gradient descent (gd) refinement. and   start-up trajectory on the strain power density chart for different admissible injector opening rate   .
.119341 Received 3 April 2023; Received in revised form 12 June 2023; Accepted 18 September 2023 [8][9][10][11][12]f the bucket (-) dynamic loads in Pelton turbines is rather restricted in the existing literature.By contrast, in reason of their wider employment, the impact of such transient operating conditions on the lifetime of Francis kg m 2 rad s ) turbines has been more deeply explored.Several publications[8][9][10][11][12]have discussed both experimental and numerical approaches assessing how transient regimes affect the fatigue of Francis turbines.Concerning  Time duration in a runner rotation during which the bucket is not hit by the jet (s)   −th time step of the start-up sequence (s)   Time duration of the jet entirely impinging the bucket in the jet charging cycle (s) . Let the polynomial objective function   (  () ,   () ,   () ) of degree   be defined by the monomials   of R[], with  ∈ N   ∶ || ≤   }, || ∶= ∑   =1   .The set of variables: [  ( 1 ) , … ,  (   ) ,   ( 1 ) , … ,   (   ) ,  1 ( 1 ) , … ,  1 (   ), ... ,  nr HF ( 1 ) , … ,  nr HF (   ) ] . The variation induced by the gradient descent refinement is lower for higher   values in reason of the moderate acceleration that characterizes the initial phase of the sequence.duration of the initial stall in the sequence with   = 2 × 10 −2 s −1 , the improvement provided to this regime is lower compared to the two others simulated sequences.Nevertheless, the solution provided by the SDP relaxed problem reaches the global optimum for all the investigated scenarios.The refinement procedure is also implemented on the shorter startups with higher injector opening rate.Figs.18and19present the improvement provided by the gradient descent method on the start-up sequences shown in Figs.9 and 10.Compared to the sequences featuring lower admissible time