Constrained optimal control of a point absorber wave energy converter with linear generator

This paper investigates a method for optimal control of a point absorbing wave energy converter by considering the constraints on motions and forces in the time domain. The problem is converted to an optimization problem with the cost function being convex quadratic and the constraints being nonlinear. The influence of the constraints on the converter is studied, and the results are compared with uncontrolled cases and established theoretical bounds. Since this method is based on the knowledge of the future sea state or the excitation force, the influence of the prediction horizon is indicated. The resulting performance of the wave energy converter under different regular waves shows that this method leads to a substantial increase in conversion efficiency.


I. INTRODUCTION
In recent years, the increasing costs of fossil fuels and the environmental problems derived from their exploitation have led to an increasing interest in producing electricity from renewable or alternative energy sources.2][3] The process of energy conversion by a WEC usually consists of several steps: wave energy absorption by the capture system, conversion of mechanical power to electricity by the power take-off (PTO) system, and transmission to the electrical grid.
Efficient conversion will be achieved when the incident wave frequency is close to the natural frequency of the WEC. 4 For the given point absorber WEC with fixed dimensions, the natural frequency is fixed, and since the sea states vary over time, optimization of the WEC is necessary.One method to do this is by automatic control of PTO parameters, which influences the performance and power production of the WEC.Damping must be adjusted to a proper level to achieve maximum energy conversion: if the damping is too high, the motions are limited and only a small amount of power is produced, but if the damping is too low, the damper does not absorb much power and most power is not captured.With any PTO system, the matching damping is vital for an efficient WEC. 1 It is well known that optimal control of the PTO in order to approach an optimum interaction between converter and incident wave can increase the power output.This has been a key issue for the researches of wave energy in the last decades.It was originated by Budal 5 and Salter 6 independently in the 1970s.Then Evans 7 and Pizer 8 studied the maximum wave power absorption under motion control.Most of the related works are performed in the frequency domain.Two widely studied approaches to control the WECs are reactive control 9 and latching control. 10Reactive control involves application of force on the device that is in phase with displacement and acceleration, and can be implemented relatively simply.However, it needs large force and large exchange of energy than what can be provided by the PTO systems usually proposed for the WECs, and cannot be applied for real time optimal control. 11Latching control requires that the velocity should be in phase with the excitation force, which can be achieved by latching the WEC when its velocity vanishes and then releasing it at a suitable time. 12o implement this in an actual system, a time domain study is necessary.A control strategy of the motion can be formulated by predicting the excitation force in the time domain. 13,14he constraints of amplitude, velocity, machinery force, and other features can be considered and included in the solution. 15A procedure to find optimal velocity and machinery force under amplitude constraints was developed and used to study how the constraints influence the optimal motion of heaving the WECs. 16Richter et al. used a state space method to present the application of nonlinear control to a point absorber WEC by considering the constraints of velocity and machinery force. 17By expanding the motions and forces in terms of truncated Fourier series, the optimal control problem of a single WEC was transformed into a constrained finite dimensional optimization problem with convex quadratic cost function. 18Later, the authors extended this application to arrays of the WECs. 19Furthermore, the problem when the model is bilinear and the cost function is not convex quadratic was solved. 20ften, control strategies for single or arrays of WECs are discussed without considering the relationship between the PTO force and the converter's motion.This leads to two problems: one is that it is difficult or impossible to implement the optimal PTO coefficients given by the numerical calculation in the true device.Another problem is that it is hard to separate the instantaneous active power from the total PTO power.In this paper, the PTO force is modelled as a function of velocity and displacement, which can be considered as another form of constraints.To simplify this, the PTO in our case consists of a spring and time-varying resistive load.Usually, the active power delivered by the resistive load is considered as the electric power produced by the converter.
We consider this as a general optimal active control problem for a heaving point absorber with a direct-driven linear generator.We follow the Fourier method introduced by Bacelli et al. 18,19 and extend it to the constrained optimal control of the WEC by considering the relationship of the PTO and the motion of the converter in the time domain.The PTO damping coefficient and the PTO spring coefficient are defined to be non-negative and possible physical constraints are included in the solution of the optimization problem.In Section II, we will discuss the hydrodynamic characteristics of a heaving WEC and formulate a model for the control strategy.Section III presents possible constraints of the WEC and studies the solution of the problem.In Section IV, the results from Section III are presented and compared with results from established theoretical upper bound and uncontrolled cases.
An established upper bound for the power-to-volume ratio, introduced by Budal and Falnes, 21 is expressed as P max =V < qgwA=4, where V is the volume of the buoy, q the density of ocean water, g the acceleration of gravity, w is the radian frequency of incident wave, and A the amplitude of incident wave.The inequality can be rewritten as P max < VqgwA=4, and P max is the maximum power produced by the WEC under this constraint, called Budal and Falnes upper bound (BFUB) in this paper.Another concept used in the comparison is the capture width (CW), which is defined as CW ¼ P mea =P avi , where P mea is the mean power consumed by the damping part of the PTO, and P avi ¼ A 2 qg 2 =ð4wÞ is the available wave power in the incident wave train per unit crest length.The maximum capture width is CW max ¼ k=ð2pÞ, which yields the maximum equivalent crest length that can be absorbed by an axisymmetric buoy in a symmetric mode of motion, e.g., heave. 5,22The relative CW (RCW), defined as the ratio of the CW and the buoy diameter, is also used to quantify the efficiency of the converter.

II. GOVERNING EQUATION
A point absorber is a device that has small dimensions compared to the incident wavelength.In this paper, we consider the WEC developed at Uppsala University which can be regarded as a point absorber with a direct-driven linear generator as PTO, as sketched in Fig. 1.The main parameters of the WEC used for calculation are shown in Table I.The linear generator is placed on the seabed and connected to a buoy via a line.In this way, the natural wave motion is transferred to the translator by the buoy's motion.For simplification, only the heave direction is considered, and the system can be modelled as a one-dimensional mass-springdamping system.
The fluid flow is assumed to be inviscid, incompressible, and irrotational.If the wave height and motion are small, the motion of the WEC can be described according to Newton's law as 24 m€ zðtÞ ¼ F e ðtÞ þ F r ðtÞ þ F h ðtÞ þ F PTO ðtÞ; ( where m is the inertia mass of buoy and translator, zðtÞ represents the vertical displacement of the buoy from equilibrium, _ zðtÞ and € zðtÞ represent the velocity and acceleration of the WEC, respectively, and F e ðtÞ is the excitation force.F h ðtÞ is the hydrostatic restoring force due to buoyancy and gravity, which is proportional to the displacement F h ðtÞ ¼ ÀSzðtÞ; (2) where S ¼ qgpr 2 is the hydrostatic stiffness coefficient with r being the radius of buoy.
The radiation force, represented by F r ðtÞ, describes the force due to the movement of the buoy itself in the absence of incident waves.One expression of it in time domain was introduced by Cummins, 24 Here, M 1 is the asymptotic value of added mass in the limit towards infinite frequency.The kernel of the convolution integral, K(t), is the radiation impulse response, and can be solved from Ogilvie's relations, FIG. 1. Schematics of the wave energy converter.Reprinted with permission from J. Appl.Phys.102, 084910 (2007). 23opyright 2007 AIP Publishing LLC.
where C(w) is the radiation damping and MðwÞ the added mass in the frequency domain.
There are other possible ways to represent the radiation force.One is introduced in Refs. 25 and 26, where the radiation force is decomposed into components in phase with the buoy's velocity and acceleration.Another method is introduced in Refs.20, 27, and 28, where the authors use the state space method to present this relationship.
For our WEC, the PTO system is a linear generator, and it can be represented by a damping force proportional to the velocity and a spring force proportional to the displacement where F PTO ðtÞ is the PTO force, K PTO is the spring coefficient, and C PTO ðtÞ is the timedependent damping coefficient that needs to be solved.Finally, Eq. ( 1) can be rewritten as

III. CONTROL AND CONSTRAINTS A. Cost function
To absorb maximum energy within the infinite time horizon, a control problem has to be formulated.This is defined as: find optimal velocity or PTO coefficients to get the maximum energy over an infinite time horizon by considering the constraints of the WEC.Since the sea states of the whole time horizon must be predicted to calculate the excitation force before the control problem, the length of the time horizon is called prediction horizon (T) here.The instantaneous power of the PTO consists of two parts: the PTO spring power and the PTO damping power.The PTO damping power is the final power produced by the converter while the PTO spring power will be returned to the ocean.In other words, the time average value of the PTO spring power is zero even though it can lead to a big flux of instantaneous PTO power.
The instantaneous power equals the product of the PTO force and the velocity.The total energy converted by the WEC over time T is Truncated zero-mean Fourier series is used to approximate the velocity and the PTO force, which leads to a finite dimension optimization problem.The basis vector is chosen as L ¼ ½cos ðw 0 tÞ; sin ðw 0 tÞ; cos ð2w 0 tÞ; sin ð2w 0 tÞ; …; cos ðNw 0 tÞ; sin ðNw 0 tÞ T ; where w 0 is the fundamental frequency of the basis and 2N is the length of the vector.It should be noted that, with this chosen control approach, the solution will depend on the choice of the basis functions, which may have influence to the optimal results.Velocity of the WEC can be expressed in terms of the basis vectors as The damping force of the PTO can be expressed as c an cos ðnw 0 tÞ þ c bn sin ðnw 0 tÞ; (11)   and similarly for the excitation force.Then an approximation of the solution to the equation of motion is solved by using the Galerkin method.The result (for more details, see Ref. 18) is where V, A, F are vectors of Fourier coefficients of velocity, PTO damping force, and excitation force, The matrix H is block diagonal and its l th block elements H l can be expressed as follows with l ¼ 1, 2… N, where By substituting Eqs.(10) and (11) into Eq.( 8), the cost function of the resulting optimization problem can be written as W ¼ A T V. 19 If the matrix H is non-singular, the optimal A that maximizes the total converted energy can be obtained by solving the optimization problem,

B. Constraints
The WEC uses a linear generator as the PTO system and the power will be produced by the relative motion of stator and translator.It is well known that the output of a system is determined by the input (such as sea states, load) and the characteristics of the system itself (such as the inertial mass and the spring coefficient).So, for a WEC without constrained control, the optimal operation of the converter sometimes may exceed the physical limitations imposed by the system (such as the amplitude of velocity, motion, or force) under certain sea states.Such situations are dangerous for the practical WEC.Especially for waves with high amplitudes and short periods, the motion of the WEC can exceed the designed stroke length, which may cause a strong strike to the hull.The improved survivability is a very interesting secondary benefit with motion control.In our device today we have motion control (limited stroke length), but the striking force can still be very high if the motion is not controlled in the correct way.

Velocity constraints
For a given buoy and a given incident wave, the excitation force is a given quantity, and the velocity and PTO force are approximated by a linear combination of the basis functions.Knowing the PTO force (velocity), the velocity (PTO force) can be found by Eq. ( 12).The constraints on the velocity and PTO force will influence each other.The constraints of velocity can be expressed as follows with v max being the maximum value of velocity and vðtÞ ¼ L T V:

Displacement constraints
Constraints on displacement can play a role in protecting the device by letting the translator move within the designed stroke length, whose value depends on the sea states where the WEC will be deployed and will be given by the designer.Two generators designed at Uppsala University are denoted L1 29 and L3 22 and these have a maximal stroke length of 1.8 m and 2.2 m, respectively.We study cases with small wave height in the present paper, and use values no larger than 0.5 m as displacement constraints.Displacement is expressed in terms of the basis vectors as With z max being the maximum value of displacement, this constraint can be written as

PTO constraints
A large and rapidly varying control force is required in a latching control, and its amplitude can be included in the solution of the optimization problem. 16The amplitude constraints on the PTO damping force can be expressed as where F PTO C;max is the maximum value of the PTO damping force and F PTO C ðtÞ ¼ L T A ¼ L T ðHA À FÞ.

Non-negative coefficients
The PTO force is split into two parts, as shown in Eq. ( 6).The PTO spring coefficient is constant in this paper.The PTO damping coefficient is the variable needs to be solved, since the spring is the coil spring obeying the Hooke's law and the damping is purely resistive in our WEC.To avoid non-physical values, the spring and damping coefficients should be defined as non-negative, which means the spring force and damping force are in the opposite direction of displacement and velocity, respectively.

Control with constraints
The issue is now converted to a minimization problem with a quadratic cost function and linear as well as nonlinear constraints.In a real-time implementation, the optimization problem is calculated over a time horizon T, called prediction horizon in this paper, with a correct estimation of the excitation force.In this case, the correct values will finally depend on the prediction of the sea states within the whole prediction horizon.Even though a long prediction horizon is desirable, it is a challenge to predict the waves over a long distance.
The constraints on the sign of the PTO damping coefficient result in a nonlinear problem, which means the issue is not a general quadratic problem.To solve this, a control code is implemented using the active-set method in MATLAB.The constraints are imposed only at specified time instants t i ¼ N i Ã dt in the range [0, T], N i are integers starting from zero and the time step is dt ¼ 0:1 s.
Hydrodynamic parameters are calculated using the boundary integral potential flow solver WAMIT.A prediction horizon of T ¼ 2p=w 0 ¼ 62:8 s and a Fourier length of N ¼ 80 are preferred in this paper, and other values are also used for comparison in Sec.IV.

Control without constraints
The optimal results without constraints, such as the velocity and mean power, can be found by solving the cost function.One example is shown in Ref. 16, where the optimal velocity is found assuming a perfect matching with and without displacement constraints using the discrete-time expression.In this paper, if the matrix H is nonsingular, the optimal damping force coefficients can also be derived by minimizing Eq. ( 17) as The optimal PTO damping coefficients in the case without constraints can then be solved for as

No control
The no control means the WEC is matched with constant loads.The displacement can be found solving the frequency-domain equation where x is the Fourier transform of zðtÞ, F is the Fourier transform of F e ðtÞ.K PTO and C PTO are both constant in this case.

A. Convergence and prediction horizon study
The results may be influenced by some parameters from the control algorithm, such as the length of Fourier series N, the fundamental frequency w 0 , and parameters from the system itself, such as the spring coefficient, the damping coefficient, the mass, and the buoy geometry.First of all, we should make sure that the algorithm is convergent.
Different values of N and w 0 are used in the calculation to check the control algorithm.The results for each case are convergent, but fluctuate with the changes of N and w 0 until both N and N*w 0 are large enough.This fluctuation due to that the motions and forces are expanded in truncated Fourier series, and certain lengths of Fourier series (N) and frequency range (N*w 0 ) are mandatory.As shown in Fig. 2, for the assumed regular wave, the produced mean power by the WEC is calculated with a different Fourier length and w 0 ¼ 0.1 rad/s or w 0 ¼ 0.05 rad/s, which corresponding to the prediction horizon of 62.8 s and 125.6 s, respectively.The results level out and the mean power is more close to the available power when the product of N and w 0 is bigger than 6.0.
With satisfying the above requirements, the influence of prediction horizon are studied as shown in Fig. 3.The amplitude of incident wave is 0.5 m and the radian frequency of incident wave is 2.0 rad/s.The full line is CW and the dashed line is the ratio of peak power P eak and mean power P mea as a function of the prediction horizon.Here, the Fourier length N and w 0 are constant, and the prediction horizon is increased from 5.0 s to 100.0 s. Results show that prediction should be long enough and at least include several wave cycles, then the CW calculated by this method will be constant.The ratio of peak power and the mean power decreases with the increasing prediction horizon until it flattens out at about 17 times the wave period, and a small prediction horizon leads to a fluctuation of the ratio.This indicates that the prediction horizon has little influence to the mean power produced by the WEC, but a small prediction horizon will lead to large power flux.

B. Spring influence
The spring coefficient in our WEC is constant, and its value will influence the performance of the converter.As shown in Fig. 4, the RCW of the WEC under different spring coefficients are calculated using the control method.Results indicate that the RCW varies with the radian frequency of the incident waves.For the same frequency, the value of the RCW will be higher if the PTO spring coefficient is smaller.In other words, a negative value of K PTO will lead to a higher RCW, which means the efficiency of the WEC will be higher and more power will be produced.This is also claimed in Refs.8 and 30.
However, a negative spring coefficient is not physically realizable in our case.In our WEC, the spring force will always be in the opposite direction to the displacement of the generator independent of its direction of motion, which requires the value of K PTO to be nonnegative in the mathematical model.Point A in Fig. 4 shows that the RCW of the WEC is about 100% when w ¼ 1.4 rad/s if negative spring coefficient is used, corresponding to the power absorbed by the WEC being equal to the available power of the incident wave.So if the FIG. 2. Results from the convergence study.The regular wave is A ¼ 0.5 m, w ¼ 1.0 rad/s.Constraints are K PTO ¼ 0; v max ¼ 1:0 m=s , z max ¼ 0:5 m, and F PTOC;max ¼ 1 Â 10 5 N.

043127-8
Wang et al.PTO is considered as an independent variable which does not have a direct relationship with the motion of the WEC, and the sign of the PTO coefficients is not considered, two problems will occur:

J. Renewable Sustainable
(1) The negative spring coefficient will be included in the optimization results.
(2) It is hard to separate the instantaneous power actually produced by the WEC from the total PTO power.
As shown in Fig. 5, the PTO power contains the PTO damping power and the PTO spring part.The damping part is always negative and its absolute value is the power produced by the FIG. 3. Influence of the prediction horizon.The y axes is non-dimensional.E ¼ 1.0 m.P eak is the peak power, P mea is the mean power, and CW is the capture width.The regular waves are A ¼ 0.5 m, w ¼ 2.0 rad/s.Constraints are  converter.The spring power is high in this case and leads to a bi-directional flow of instantaneous PTO power, which means that the PTO returns some power to ocean.

C. Influence of constraints and control
Many constraints can be imposed on the WEC.These constraints, such as velocity, displacement, and the PTO force constraints may be closely connected, and influence the performance of the WEC together.It is hard to study the factors separately, so we use several cases with different velocity and displacement constraints under regular waves with wave amplitude of 0.5 m to evaluate the influence of the constraints and performance of the WEC.The mean power produced by the WEC under such cases is compared with the available power in the wave.
As shown in Fig. 6, when the displacement constraint is kept constant, and the velocity constraint is increased from 0.5 m/s to 1.0 m/s and then to infinity, the WEC exhibits little difference in short-period waves and no difference in long-period waves.However, the mean power will be highly decreased by changing the displacement constraint from 0.5 m to 0.2 m while keeping the velocity constraint constant, or by changing the velocity constraint to a very low value (0.1 m/s).It seems that both the velocity constraint and displacement constraint have a great influence on the result.This is due to the fact that this control method requires the displacement of the WEC to increase quickly to a proper value at some moments, which is similar to the latching control.This behavior will be destroyed if the velocity or displacement amplitude is limited to small values, resulting in less captured power.Furthermore, the WEC has good performance when v max ¼ 1:0 m=s and z max ¼ 0:5 m, better than other five cases, and almost reaches the same mean power as in the case where the velocity constraint is infinity.Also, considering the stroke length of our practical WEC, it indicates that v max ¼ 1:0 m=s and z max ¼ 0:5 m are realistic constraints on our WEC.Fig. 7 shows the displacements of the WEC under three controls, namely, control with constraints, control without constraints and no control.One obvious phenomenon is that the phase of the displacement under control lags the uncontrolled case.The displacement amplitudes are also quite different: the displacement amplitude from uncontrolled case is lower than the incident wave height (0.6 m), while the optimal displacement amplitude for the case without constraints is more than four times the wave height.

043127-10
It is well known that the WEC will run with optimal motion and produce more power with proper load matching to the expected sea state.The results in Fig. 8 show that the WEC has a good performance and can produce at least 2.0 kW power across a large frequency range of 0.2-2.0rad/s using the control method, which is better than the uncontrolled cases.The WEC under no control has a narrow bandwidth, and its performance varies with the loads.9 shows the motions and forces of the WEC with considering the constraints.The constraints are v max ¼ 1:0 m=s and z max ¼ 0:5 m, and the wave height is 0.6 m.It indicates that the excitation force is in phase with the velocity, and the displacement seems to be latched when the velocity vanishes, which is more obvious for the long-period wave, as shown in Fig. 10.In general, the control method leads to a performance that bears resemblance to performance under latching control and produces more power.

D. Performance of the WEC in different regular waves
By using the constrained optimal control method introduced in this paper, the performance of the WEC in different regular waves with wave amplitude of 0.3 m, 0.4 m, and 0.5 m is evaluated in this section.The mean power produced by the WEC with constraints of v max ¼ 1:0 m=s and z max ¼ 0:5 m is compared with the BFUB in Fig. 11.It indicates that the produced mean power is close to or slightly exceeds the BFUB at low frequency.It should be noted that different displacement and velocity constraints are used in the BFUB, which has great influence on the mean power as introduced in Sec.IV C. In their derivation of the upper bound, Budal and Falnes assumed that the displacement amplitude does not exceed V/2S w and the velocity amplitude does not exceed Vw/2S w , where S w is the water plane area. 21Here, this assumption does not hold due to the non-harmonic velocity of the buoy, and constraints realistic for our WEC are used.As pointed out previously, the chosen control approach or the basis function may also influence the optimal results.Pizer 8 and Hals 31 also reported that the power absorption can exceed the limit bound under the additional restriction of sinusoidal motion.
The corresponding capture widths are shown in Fig. 12.The capture width never exceeds the maximal capture width, and it is close to the buoy's physical width when w ¼ 1.4 rad/s, which means the WEC has a high RCW exceeding 90%.The ratio of peak power and mean power are used to evaluate the peak variation of power in different regular waves.Results show that the ratio decreases with the radian frequency for the constant prediction horizon.Since the radian frequency is inversely proportional to the wave period and the prediction horizon is constant, it means that the ratio will decrease if more wave cycles are considered in the constant prediction horizon.In other words, a long prediction horizon can decrease the variation of power.Another case is studied as shown in Fig. 3 with changing the prediction horizon while keeping the wave amplitude and wave height as constant.Results show that the ratio decrease with the prediction horizon, and peak power will be several times of the mean power in long prediction horizon.
However, it is hard to have an infinite prediction horizon in the practice.Since this method is based on the model predictive control, which means the future excitation force used in the prediction horizon must be predicted exactly by measuring waves at a certain distance before the WEC.So it is necessary to have proper prediction horizon for the optimal control of the WEC.

V. CONCLUSIONS
A time-domain analysis has been presented to evaluate the performance of a point absorber WEC whose motion is governed by a control algorithm.The WEC is oscillating in heave direction and in regular waves.By considering constraints of motion and forces, a convergent numerical scheme was implemented to solve the power maximization problem.It indicates that constraints can influence the performance of the converter: with suitable constraints, less losses, and more absorption power can be achieved.In the method, future knowledge of the excitation force is assumed within a relatively short finite time interval (the prediction horizon).
It was found that, although no phase control has been directly attempted, in theory this control method can attain a phase shift between excitation force and velocity of buoy, which means the excitation force is in phase with the velocity, and this performance is similar to the FIG.12. Performance of the WEC in different regular waves (amplitudes of incident waves are 0.3 m, 0.4 m, 0.5 m, respectively).P eak /P mea means the ratio of peak power and mean power.CW is the capture width, CW max is the maximum capture width, and they are non-dimensioned by a length scale E ¼ 1.0 m. performance from latching control.High relative capture width, more than 90% in some regular waves, can be achieved by this method.

FIG. 4 .
FIG. 4. Relative capture width (RCW) of the WEC with different PTO springs.The wave amplitude of the incident wave is 0.3 m. z max ¼ 0.5 m, v max ¼ 1.0 m/s.

FIG. 6 .
FIG.6.Performance of WEC with different constraints.The amplitude of incident wave is 0.5 m.Fundamental frequency is 0.1 rad/s.PTO spring coefficient is zero.

Fig.
Fig.9shows the motions and forces of the WEC with considering the constraints.The constraints are v max ¼ 1:0 m=s and z max ¼ 0:5 m, and the wave height is 0.6 m.It indicates that the excitation force is in phase with the velocity, and the displacement seems to be latched when the velocity vanishes, which is more obvious for the long-period wave, as shown in Fig.10.In general, the control method leads to a performance that bears resemblance to performance under latching control and produces more power.

FIG. 8 .
FIG. 8. Mean power of the WEC with and without control.The wave amplitude is 0.3 m, and z max ¼ 0.5 m, v max ¼ 1.0 m/s, K PTO ¼ 0.

FIG. 10 .
FIG. 10.Motions and forces of the WEC under control with constraints.The regular wave is A ¼ 0.3 m, w ¼ 0.4 rad/s.Constraints are z max ¼ 0.5 m, v max ¼ 1.0 m/s, K PTO ¼ 0.

TABLE I .
Main parameters of the WEC.