Model Predictive Control of CorPower Ocean Wave Energy Converter

(1)

INOM TEKNIKOMRÅDET EXAMENSARBETE

TEKNISK FYSIK

OCH HUVUDOMRÅDET ELEKTROTEKNIK,

AVANCERAD NIVÅ, 30 HP STOCKHOLM SVERIGE 2016,

Model Predictive Control of CorPower Ocean Wave Energy Converter

CAROLINA ERIKSSON

(2)

Model Predictive Control of CorPower Ocean Wave Energy Converter (WEC)

Carolina Eriksson careriks@kth.se

16 juni 2016

(3)

Abstract

Wave power is currently a hot topic of research, and has shown great potential as a renewable energy source. There have been lot of progress made in developing cost effective Wave Energy Converters (WECs) that can compete with other sources of energy in regard to price and electrical power. Theoretical studies has shown that optimal control can increase the generated power for idealized WECs. This thesis is done in collaboration with CorPower Ocean, and investigates the use of economic Model Predictive Control (MPC) to control the generator torque in a light, point-absorbing, heaving WEC that is currently under development. The objective is to optimize the generator torque, such that the average generated power is maximized while maintaining a small ratio between maximum and average generated power. This results in a nonconvex cost function. Due to the highly nonlinear and nonsmooth dynamics of the WEC, two controllers are proposed. The first controller consists of a system of linear MPCs, and the second controller is a nonlinear MPC. Relevant forces acting on the WEC are identified and the system dynamics are modelled from a force perspective. The models are discretized and the controllers are implemented in Simulink. The WEC, together with the controllers, is simulated in an extensive Simulink model developed by CorPower Ocean. Several different types of ocean waves are considered, such as its energy content and its regularity.

In the majority of cases, the controllers do not increase the performance of the WEC compared to a simple, well tuned controller previously developed by CorPower Ocean. Finally, possible improvements of how to reduce existing model errors are proposed.

(4)

Sammanfattning

V˚agkraft har de senaste ˚aren visat stor potential som en ny, förnyelsebar energikälla. Det har skett m˚anga framsteg inom omr˚adet med att ta fram ett robust v˚agkraftsverk som kan utmana andra energikällor i pris och elektrisk effekt. Teoretiska studier har visat att optimal styrning kan öka den elektriska effekten för idialiserade v˚agkraftsverk. Denna rapport är skriven i sammarbete med v˚agkraftföretaget CorPower Ocean, och undersöker hur ekonomisk Model Predictive Control (MPC) kan användas för att styra dämpningen i ett lätt v˚agkraftverk vars storlek är relativt liten v˚aglängden. M˚alet är att optimera dämpningen, vridmomentet, i generatorerna s˚a att medeleffekten maximeras samtidigt som toppeffekten minimeras, detta för att skapa ett stabilare system med mindre flutuationer mellan medel- och toppeffekt. För att n˚a detta m˚al krävs en icke konvex kostfunktion. P˚a grund av stora olinjäriteter och diskontinuteter i systemets dynamik utvecklas tv˚a regulatorer; ett system av linjära MPC, samt en olijär MPC. Relevanta krafter som p˚averkar systemet identifieras och modelleras fr˚an ett kraftperspektiv. Modellerna diskretiseras, och regulatorerna implementeras och simuleras i en detaljerad Simulink modell av systemet, utvecklad av CorPower Ocean. B˚ade regelbundna och oregelbunda v˚agset med varierande energiinneh˚all har simuleras. Regulatorerna ökar inte v˚agkraftverkets prestanda jämfört med en enkel, väl inställd regulator utveklad av CorPower Ocean. Slutligen föresl˚as förbättringar för att minska modelfell i modellerna.

(5)

Acknowledgements

I want to express my gratitude towards CorPower Ocean for letting me write my thesis together with them, and letting me be part of their exciting project. Thanks to my supervisors Gunnar Steinn ´Asgeirsson and Jørgen Hals Todalshaug, for their enthusiasm, assistance and for giving insightful inputs whenever I needed it. Thanks to Elisabet Jansson, for answering my all my questions regarding the Simulink model, both small and large. Special thanks to my examiner at KTH, Elling Jacobsen for his valuable encouragement and advise. Finally, I would like to thank my partner Hjalmar Heimb¨urger, for his unremitting support, and patience during this project.

(6)

List of Figures

2.1 Coordinate system of the WEC.. . . . 4 2.2 The point absorber interfere with the incoming waves by oscillat-

ing in the water. The radiated waves depends on the amplitude and frequency of the buoy [1]. . . . . 5 3.1 Forces acting on the WEC. . . . . 13 3.2 The ideal case of flywheel engagement and disengagement to the

rack of the WEC.. . . . 15 3.3 Efficency η of the generators, dependent on the velocity of the

flywheel and the torque in the generator. Data provided by CorPower Ocean. . . . . 16 3.4 The net force Fnet= Fhydrostatic+ Fwave spring+ Fpretensionfrom

hydrodynamic and phase control. Data provided by CorPower Ocean. . . . . 17 4.1 The connection between the engaged LMPC and the disengaged

LMPCs. . . . . 19 4.2 F_netwith two different estimations of ˆF_net, y_linear and y_poly. . . 20 4.3 Estimation of the efficiency. The point close to zero have been

added to ensure that the function passes close to origo.. . . . 21 6.1 Normalized mean power output for regular waves, sea state 1,

LMPC. . . . . 32 6.2 Normalized mean power output for regular waves, sea state 2,

LMPC. . . . . 33 6.3 Normalized mean power output for regular waves, sea state 3,

LMPC. . . . . 34 7.1 Above: Simulated position and velocity of the rack of the WEC,

for sea state 1. Below: The observed behavior of the switchnig for sea state 1 (blue) compared to the expected switching behavior (red). . . . . 39 7.2 Expected vs. real position and velocity of WEC, Model errors

causes a large decrease in velocity instead of a small increase. . . 41 7.3 Expected vs. real position and velocity of WEC. Weird behavior

for velocities close to zero. . . . . 42

(9)

List of Tables

5.1 Sea State Parameters. . . . 28

6.1 Reference values for each sea state, regular waves. . . . . 29

6.2 Reference values for each sea state, irregular waves.. . . . 29

6.3 Linear Damping, sea state 1, with maximum damping force v_max= 5400 · k_gen N and h = 0.2, regular waves. . . . . 30

6.4 Mean power [kW] for regular waves, sea state 1, LMPC. . . . . . 30

6.5 Ratio of maximum to mean power, for regular waves, sea state 1, LMPC. . . . . 31

6.10 Mean power [kW] for irregular waves, sea state 1, LMPC. . . . . 35

6.11 Ratio of maximum to mean power, for irregular waves, sea state 1, LMPC. . . . . 35

6.16 NMPC simulated with regular waves for sea state 1. . . . . 37

6.19 NMPC simulated with irregular waves for three different sea states and different sample time, N = 5. . . . . 37

7.1 Mean Power output [kW], LMPC, with perfect switching for regular waves. Sea state 1.. . . . 40

7.2 Peak to Average Ratio, LMPC, with perfect switching for regular waves . . . . 40

Nomenclature

(10)

A State space matrix

B State space matrix

C State space matrix

Frad tot The total radiation force Fhydrostatic The hydrostatic force Fdrag The drag force

Fwave spring The force from the WaveSpring Fpretension The force from the Pretension cylinder

Fnet The sum of the WaveSpring, Pretension and hydrostatic forces FP T O The force from the Power Take-Off system

Fhydro The total hydrodynamical forces Fmech The total mechanical forces F_exc The excitation force f₀(t, x, u) The cost function

f (t, x, u) General system dynamics

g Gravity

H_s Significant wave height

h Sample time

i Index of the generators and flywheels, i ∈ [1, 2]]

J Inertia

kgen Conversion ration between translational to rotational movement N Number of prediction steps into the future

η Efficiency

Φ(t) Terminal cost

q Radiation state vector ρ Mass density of water Tp The peak period

t Time

u(t) General control signal V Set of allowed control signals V_sub Submerged volume of WEC v_i Torque for the i:th generator

v(kh) Discretized zero order hold control signal wexc Excitation force

X Set of allowed state values x(t) Dynamic state vector xd(k) Discrete state vector ξ˙_i Velocity of flywheel i ξ¨_i Acceleration of flywheel i y General output vector

z Buoy position, in the heave direction

˙

z Buoy speed, in the heave direction

¨

z Buoy acceleration, in the heave direction

(11)

Abbreviations

CPO CorPower Ocean

WEC Wave Energy Converter MPC Model Predictive Control PTO system Power Take-Off system LTI Linear Time Invariant

(12)

Chapter 1

Introduction

1.1 Background

There is a need to replace fossil fuels with more sustainable alternatives, and reduce the world’s CO₂ emissions. Sustainable alternatives such as solar and wind power are well developed technologies that are used around the world.

Wave power is not as developed as solar and wind technologies [2] and it is calculated that the oceans can deliver around 1 – 10 TW of energy [3]. The oceans is therefore an important unused source of sustainable energy and an interesting research area.

There exist several different proposals on potential Wave Energy Converters, WECs, such as tidal energy converters, oscillating water columns and surface heaving buoys. This thesis will focus on a light, point-absorber heaving buoy, that is currently being developed by Corpower Ocean AB. The buoy absorbs energy from the surface waves by resonating with the waves, while being moored to the sea bed.

Previous research has shown that the use of phase control can make WECs oscillate in resonance with the incoming waves to create the most optimal velocity pattern [1]. Several different ways of controlling the phase of the buoy have been proposed, such as passive control (adapting the damping to the current seastate) and latching [4]. Corpower Ocean have been developing the innovative WaveSpring, a system with negative springs that creates resonance with the incoming wave and thus convert more wave energy into electrical energy [5].

When an optimal phase control has been achieved, the generator torque needs to be controlled. This is because the generated power depend on the generator torque and on the rotational speed of the generator. The generators act as dampers that reduce the movement of the buoy, and the choice of torque can therefore also affect both the speed of the generators and the velocity of the WEC. It is therefore important to find the most optimal choice of torque in the generators such that the power produced by the WEC can be maximized.

(13)

Optimal control and model predictive control (MPC) of the generator torque have been investigated earlier for generalized and linearized WECs [6, 7, 8].

These studies found that MPC improves the performance of the WEC and increases the average power produced. However, one problem with MPC and receding horizon optimization is that the real systems rarely are linear, and that the computational time needed for solving the optimization problem is too long. There have been some research on nonlinear MPC for WECs. For example investigated Ritcher[8] the performance of a non linear MPC for model of a WEC with a highly non linear mooring force and compared it with a linear MPC. However, Ritcher did not consider the computational time, and the MPC formulation will differ depending on how the WEC is designed.

This M.Sc. thesis develops one linear and one non linear MPC for a lightweight point absorber WEC, with focus on modelling the nonlinear and nonsmooth in the Power Take-Off (PTO) system in the WEC, specially designed by CorPower Ocean.

1.2 CorPower Ocean AB

CorPower Ocean is a Stockholm based company that currently is developing a lightweight and WEC. The innovative system is believed to produce five times more energy compared to other WECs, within the same size range. The CorPower Ocean WEC consists of a PTO system and the specially designed WaveSpring, both embedded into a buoy with a pretension mooring line connecting the buoy with the sea bead.

1.3 Research Objectives

The research objectives of this thesis is to investigate the use of a predictive controller for deciding the optimal generator torque. The goal is to maximize the average generated electrical power and minimize the ratio between the peak and average generated electrical power. The controllers will be simulated in Simulink and compared with a simple reference controller designed by CorPower Ocean.

1.4 Thesis Overview

First, the theoretical background to wave power and MPC is given in chapter 2. Secondly, the CorPower Ocean WEC is described in chapter 3, with focus on the forces influencing the movement of the buoy and the dynamic of the PTO system. Two different controllers are designed for the system in chapter 4, with one linear MPC and one nonlinear MPC. The implementation and simulation of the controllers are briefly described in chapter 5, before presenting the simulation results in chapter 6. Chapter 7 discusses possible improvements

(14)

of the system, and chapter 8 summaries the report and propose possible further works.

(15)

Chapter 2

Theoretical Background

This chapter will give an theoretical background to wave power and control theory.

2.1 Coordinate system

The Wave Energy Converter, WEC, moves in six degrees of freedom, as shown inFigure 2.1. Following the notation in most works regarding wave energy, the defined coordinate system is, in following order; surge, sway, heave, roll, pitch and yaw. Surge (the x-axis) is the direction of the wave.

Only the heave movement of the buoy, the movement along the z-axis, will be considered when modelling the controllers. In this thesis z is the position of the buoy above equilibrium, ˙z is the velocity in the heave direction of the buoy and ¨z is the heave acceleration. The controllers will be simulated in a Simulink model with two degrees of freedom, surge (x-axis) and heave (z-axis).

The controllers have been implemented such that it is possible to simulate them

surge (x)

sway, (y) heave, (z)

yaw

pitch roll

Figure 2.1: Coordinate system of the WEC.

(16)

(a) Radiated circular waves interfering with incoming waves

(b) a) normal waves, b) & c) an point absorber heaving/surging therefore creating waves, d) the sum of a), b) & c) causing destructive interference, total absorption.

Figure 2.2: The point absorber interfere with the incoming waves by oscillating in the water. The radiated waves depends on the amplitude and frequency of the buoy [1].

within a Simulink model with six degrees of freedom but this is not within the scope of this thesis project and will therefore not be considered.

2.2 To Create Energy from Waves

The energy in ocean waves comes from the solar energy, that have created winds that blow over the ocean and create waves. The wave energy can travel long distances, and the energy density in waves vary with season, with up to six times more energy during the winter than during the summer [3].

To make a general WEC absorb energy from the waves means that energy need to be removed from the wave. Hence, there need to be a reduction of the waves passing by or being reflected by the WEC. By making the WEC oscillate, it can radiate waves that cancel the parts of the incoming waves, seeFigure 2.2a. The WEC will thereby absorb that part of the energy that were previously contained by the wave. This oscillation and radiation of waves means that “a good wave absorber must be a good wave maker” [1].

For optimal absorption, destructive interference between the WEC and the incoming wave is needed. To create optimal destructive interference, the buoy need to oscillate with correct amplitude and phase. Two objects of different size can create waves with the same size, that is amplitude and period, if the smaller object oscillates with a larger amplitude [1]. This property is used within wave conversion, with a small floating point absorber WEC oscillating with a higher amplitude than the incoming wave. Figure 2.2billustrate an idealized absorption rate of 100%.

(17)

2.2.1 Wave Spectrums

Deep ocean waves can be modelled and described using either regular waves or irregular waves. Regular waves is described by one single sinusoidal function with a fixed amplitude and frequency. This model does not capture the randomness and irregularity that exsist in real ocean waves. Irregular waves are modelled as a superposition of multiple sinusoidal components by using a Fourier series. The surface elevation can be expressed as [9]

h_{surf ace}(t) =

N

X

k=1

A_kcos(ω_kt + _k). (2.1)

Here is t the time, _k are random phases between 0 and 2π and A_k are Rayleigh distributed amplitudes. These are defined using the mean squared amplitude from the wave spectrum S(ω) [9]:

E[A²_k] = 2S(ω_k)∆ω_k. (2.2) With ∆ωk , ω^k − ωk−1. In general we can assume that the sea surface is stationary for time periods between 20 minutes to 3 to 6 hours. These sea sea surfaces, and corresponding wave spectrum, can be described using a set of parameters, such as the significant wave height Hs and the peak period Tp. We will in this thesis only use spectrums with only one peak period. The significant wave height is defined as the average height from trough to crest, for the highest one third of the waves in the selected time period. The peak period is the wave length corresponding to the frequency at which the wave energy spectrum got its maximum value [9]. Here, the JONSWAP wave spectrum will be used to define S(ω), together with the sea state shaping parameters H_s and T_p. The JONSWAP spectrum is a commonly used wave spectrum, and a more detailed explanation of this spectrum can be found in [9].

2.3 Control Theory

2.3.1 Sampling of Dynamical Continuous Systems

A general continuous system can be modelled as a continuous dynamic system.

The controller will act as a discrete system, and therefore must the continuous model be transformed into a discrete model. This can be done using zero order hold, or numerical methods such as forward Euler. Different types of systems are discretized with different methods. Here, both linear time invariant (LTI) systems and more general nonlinear systems are treated.

Linear Time Invariant Systems

The simplest dynamical system is a linear time invariant (LTI) system, and they are therefore easy to use when creating simple models. A LTI system can be

(18)

described as [10]

˙

x(t) = Ax(t) + Bu(t) (2.3)

y(t) = Cx(t) + Du(t), (2.4)

Where x(t) ∈ Rⁿ is the state vector, u(t) ∈ R^m is the control signal and y(t) ∈ R^p the output vector. A, B, C and D are constant matrices with dimensions n × n, n × m, p × n and p × m. This ensures valid matrix operations producing the correct output. The solution to (2.3) can be expressed as [10]

x(t) = e^A(t−t⁰⁾x(t0) +

t

Z

t₀

e^A(t−s)Bu(s)ds. (2.5)

To sample the system to a discrete system, let h be the sampling time.

Integrating (2.3) from the time kh to kh + h, k ∈ N, gives

x(kh + h) = e^Ahx(kh) +

kh+h

Z

kh

e^A(kh+h−s)Bu(s)ds. (2.6)

Applying zero order hold [11] and thereby restricting the control signal u(t) to be piecewise constant as

u(t) = v(k), for t ∈ [kh, kh + h) (2.7) gives a simplification to the integral above

kh+h

Z

kh

e^A(kh+h−s)Bu(s)ds =

h

Z

0

e^A(s)dsBv(k). (2.8)

Let x_d(k) , x(kh). (2.6) can now be formulated as results in the linear time invariant discrete system

xd(k + 1) = Adxd(k) + Bdv(k), (2.9)

where Ad = e^Ah and Bd =

h

R

0

e^A(s)dsB. The discrete time system obtain the same values as original system at the sampling points t = kh.

Forward Euler

For a more complex system than a LTI system, such as a general nonlinear system, is often impossible to find an analytic solution to the differential equation that can be discretized. Numerical methods can instead be used to find a numerical solution and a discretization of the system. Examples of numerical methods include forward Euler, the midpoint method, Runge-Kutta methods, or linear multistep methods. When choosing an discretization method the stability of the system is highly important. In this thesis Forward Euler will be used.

(19)

Forward Euler is a simple explicit method, that can be derived from a simple Taylor expansion of the system [12], transforming an general nonlinear system

˙

x = f (t, x(t), u(t)) (2.10) subj. to: x(0) = x₀, x ∈ X, u(t) ∈ U (2.11) (2.12) into the discrete representation, with h is the sample time,

x(kh + h) = x(kh) + hf (x(kh), u(kh)). (2.13) The biggest advantages of using Forward Euler is that it is very easy to solve and demand little computational power. However, if a too big step length h is chosen the solution will diverge from the exact solution. With a small h, the solver also need more computational power.

2.3.2 Model Predictive Control

Model predictive control (MPC) is an advanced control method that optimizes the current time step subject to the predicted future. This is achieved by at every time step solving a finite time constrained optimization problem and only applying the first optimal control signal. The algorithm can be described as follows [13, p. 423]:

1. At time k, measure, or estimate, the state vector.

2. Calculate the optimal control signal for the N next steps.

3. Apply the first control signal v(k).

4. Wait for the next time step k + 1 and go to 1.

The MPC algorithm incorporates feedback that reduces the impact of model errors and measurement noise, compared to open optimal control. It is also possible to predict future actions or changes in the system. This makes it possible to act before the changes happen, which is not possible if for example a PID regulator is used. The general MPC algorithm can be used on almost any dynamic system, whether it is linear or nonlinear. Hence, there are several different ways of implementing the controller.

In general, state estimation is necessary in step 1 in the algorithm for estimating states that cannot be measured, and to reduce measurement noise of the measured states. This is often done by using a Kalman filter [10].

(20)

Linear Model Predictive Control with Quadratic Cost Function The optimal control problem for a discrete LTI system with a finite horizon can be expressed as

min

v(k)∈V

Φ(x_d(N )) +

N −1

X

k=0

hf₀(k, x_d(k), v(k)) (2.14)

subj. to







xd(k + 1) = Adx(k) + Bdv(k) x_d(0) given, x_d(k) ∈ X v(k) ∈ V

. (2.15)

f0(k, xd(k), v(k)) is the cost function and Φ(xd(N )) is the terminal cost. The cost function and with the terminal cost defines, together with the constraints, the optimal solution v^∗(k) for k = [0, N − 1]. The optimal solution can for example correspond to minimizing deviation from a predefined state trajectory or minimizing the control signal, given constraints on x(k) [14]. One commonly used control problem is the finite discrete Quadratic Problem (QP), with a LTI dynamic and where the cost function is defined as

f0(k, xd(k, v(k))) = x^T_d(k)Qxd(k) + v(k)^TRv(k), (2.16) where Q is symmetric and positive semidefinite and R is symmetric and positive definite. The terminal cost for the QP is

Φ(x_d(N )) = x_d(N )^TQ₀x_d(N ), (2.17) Q₀being symmetric and positive semidefinite. Finally, the constraints on xd(k) and v(k) can be written as

Fxd(k) v(k)

≤ b. (2.18)

The terminal cost and terminal constraints can be used to ensure stability for the system. A converging system would otherwise not always be guaranteed.

This is a convex optimization problem for which it is possible to solve analytically for smaller systems, see [14, pp. 22-24], or for larger problems numerically by using for example fmincon in MATLAB.

Nonconvex and Nonlinear Model Predictive Control

Nonlinear Model Predictive Control (NMPC), is not as researched as the more simple linear MPC. However, the research area has grown in the past years.

Sometimes it is more valuable to directly optimize the performance of the process instead of optimizing how the process follows a given trajectory. This is known as Economic MPC, and can have a nonconvex and nonlinear cost functions [15]. Economic MPC have been proven useful in several application areas such as optimizing cooling systems for commercial refrigeration [16], or for collision avoidance for networked vehicles [17].

(21)

The general optimal control problem can be formulated as min

v(k)∈V

Φ(x_d(N )) +

N −1

P

k=0

f₀(k, x_d(k), v(k)) (2.19)

subj. to







x(kh + h) = x(kh) + hf (x(kh), v(kh)) x(0) = x0 x(k) ∈ X

v(k) ∈ V

(2.20)

with f0(·) being the economical function to optimize, for example the produced power, and Φ being the terminal constraints. The terminal constraints guarantee that the solution converge to the wanted value. With a nonconvex cost function even small problems needs to be solved numerically and it is not guaranteed that the obtained optimum corresponds to the global optimum. The solution depends very much on the solver used and on the starting guess supplied to the numerical solver. The solver needs to be faster than the chosen step length, for the controller to be used in real time situations. This is one major limitation of NMPC, since the controller is useless if the time needed for solving the optimization problem is larger than the time until the next iteration.

(22)

Chapter 3

The CorPower Ocean Wave Energy Converter

This chapter will present the general design of the CorPower WEC, from a force perspective. The hydrodynamic and mechanical forces acting on the buoy will be presented. The focus will be on the forces needed for for the modelling of the dynamic system.

3.1 Hydrodynamic Forces

The movement and oscillation of the buoy is affected by hydrodynamic forces.

These forces consists of an excitation force from incoming waves, a radiation force due to radiate waves, a hydrostatic force and an drag force due to hydrodynamic effects.

3.1.1 Excitation and Radiation Forces

An object moving in the water due to incident waves will create radiated waves, as discussed in section 2.2. These radiated waves will in turn affect the moving object. These forces that affect the object and that are created by the movement are called radiation forces and are illustrated inFigure 2.2a.

In 1962, Cummins [18] showed that the radiation forces in the time domain can be described as an convolution between the impulse-response of the radiation force and the velocity of the moving object,

F_{rad tot}= m_∞¨z +

t

Z

0

K(t − s) ˙z(s)ds. (3.1)

The infinitive-frequency added mass m_∞ represents an instantaneous force response of the fluid after an impulsive movement of the buoy. The impulse- response K(·) depend on the the buoy and can be computed using time-domain

(23)

simulations with for example the software WAMIT. Computing the convolution integral is time consuming and the convolution is therefore often estimated as a state space representation,

˙

q = Arq + Brz˙ (3.2)

F_rad= C_rq + D_rz,˙ (3.3)

where q is the radiation state space vector and F_{rad tot}can now be expressed as

F_{rad tot}= m_∞z + F¨ _rad. (3.4)

3.1.2 Hydrostatic and Drag Forces

The hydrostatic and drag forces are due to fluid dynamic effects. The hydrostatic force, also called the gravitational restoring force, strive to press the buoy back to its equilibrium. The force is proportional to the submerged volume V_sub of the buoy, and can in the heave direction be expressed as

Fhydrostatic= (ρV_sub− m_b)g, (3.5) where ρ is the density of the water and g the gravity constant. The drag force is a quadratic force that depend on the relative speed between the buoy and the water [19] and is given by

Fdrag= 1

2ρAsubCdz˙². (3.6) Here Asub is the cross sectional submerged area of the buoy perpendicular to z, and Cd is an damping coefficient that is found experimentally.

3.2 Mechanical Forces of the Corpower Ocean WEC

The Corpower Ocean WEC consists of a power take off system (PTO system) for transforming the up and down movement of the buoy into electrical power, and a hydraulic system used to create phase control and enhance the movement of the buoy. These two system create mechanical forces that will affect the movement and dynamic of the buoy. A schematic figure of these mechanical forces is showed inFigure 3.1.

(24)

Figure 3.1: Forces acting on the WEC.

3.2.1 Power Take Off (PTO) - System

The power take off system is the mechanical system that transform the vertical movement to electric energy. It consists of two flywheels connected to the rack through a gearbox and two generators, each connected to a flywheel.

The flywheels act as energy storage and switch between being connected (engaged) to the rack, and being disconnected (disengaged) to the rack. This switch is done by a freewheel that is placed between the flywheel and the gearbox. When the flywheel is engaged to the rack, the rack and the flywheel moves with the same acceleration as the rack, and the same relative rotational speed. When the rack start to decelerate the flywheel gets a higher speed than the driveshaft from the gearbox, and disconnect. When disconnected, the

(25)

flywheel is only damped by the inertia of the flywheel, the friction in the system, the generator and by the torque applied to the generator. If the friction in the system is neglected, the dynamic for the flywheel can be described as,

ξ¨i=

(−_J ^vⁱ

f w+Jgen if disengaged

kgenz¨ if engaged , (3.7)

where ¨ξi is the acceleration of the i:th flywheel and vi is the applied torque on the i:th generator, with i ∈ [1, 2]. The constant kgenis conversion ratio between translational to rotational movement, given by the size of the cogwheels in the gear box. The ideal switching take place when the flywheel switch from engaged to disengaged when the rack reaches extreme velocities (maximum/minimum), and engage when the velocity of the rack switches sign again three quarters of the oscillation period later. This is illustrated inFigure 3.2.

The damping force acting on the buoy from the PTO system can be modelled as:

FP T O=

(−Jgear boxzk¨ gen if no flywheel is engaged

− (−sign( ˙z)(Jf w+ Jgen) + Jgear box) kgenz − k¨ genvi if flywheel i is engaged (3.8) When no flywheel is engaged, the buoy is only damped by the inertia of the

gear box. When a flywheel is engaged the buoy is damped by both the inertias from the generator, the flywheel and the gear box, as well as the torque applied in the generator.

(26)

0 0.2 0.4 0.6 0.8 1 1.2 1.4

−1.5

−1

−0.5 0 0.5 1 1.5

Time [t/T_p]

Normalizedresponce z [m]

˙ z [m/s]

Flywheel 1 disengaged Flywheel 2 disengaged Flywheel 1 engaged Flywheel 2 engaged

Figure 3.2: The ideal case of flywheel engagement and disengagement to the rack of the WEC.

The flywheels are in turn connected to generators that transform the rotational energy of the flywheels to electrical power. The damping of the generator can be controlled by changing the generator torque. The generators are not ideal and therefore there will be losses in the generator due to friction. The generator efficiency vary with the torque and the rotational speed, and the electrical power output [20] can be calculated as

P =

T

Z

0

v1ξ˙1· η( ˙ξ1, v1) + v2ξ˙2· η( ˙ξ2, v2)dt, (3.9)

where η denotes the efficiency of the generators. The efficiency of the generator is visualized inFigure 3.3.

(27)

Figure 3.3: Efficency η of the generators, dependent on the velocity of the flywheel and the torque in the generator. Data provided by CorPower Ocean.

3.2.2 Phase Control

The Corpower WEC have designed the inovative WaveSpring to achieve phase control of the buoy. The WaveSpring consist of two negative springs that have been shown to reduce the hydrostatic force [5]. The WaveSpring can be tuned to achieve optimal phase control, together with a pretension spring. 3.4show the total net force when adding up the force from the WaveSpring, the pretension spring and the hydrostatic force.

(28)

−6 −5 −4 −3 −2 −1 0 1 2 3 4 5 6

−4,000

−2,000 0 2,000 4,000

z [m]

Force[kN]

F_WaveSpring Fhyst

Fgas

FWaveSpring+ Fhyst+ Fgas

Figure 3.4: The net force F_net = Fhydrostatic+ Fwave spring+ F_pretension from hydrodynamic and phase control. Data provided by CorPower Ocean.

3.3 Kinematic Relationships

We have divided the forces that act on the buoy into hydrodynamic forces (section 3.1) and mechanical forces (section 3.2). Together with Newton’s second law these forces can be combined into the differential equation

M ¨x = Fmech(t) + Fhydro(t) (3.10) where x is the state vector and M the inertia of the buoy. The forces acting on the system is shown in Figure 3.1, and can be summarized into the following equations

F_hydro= F_exc+ F_drag+ Fhydrostatic+ F_rad (3.11) F_mech= F_trans+ F_pretension+ F_{f riction}+ F_wavespring. (3.12) By replacing Fhydrostatic, F_pretension and F_wavespring with F_netfromFigure 3.4 and expressing the excitation force as

F_exc= w(t), (3.13)

(3.14) can be expanded and the kinematic relationships of the WEC can be written as

M ¨x = F_trans+ F_{f riction}+ F_net+ F_drag+ F_rad+ w(t). (3.14) With the transmission force from (3.8), excitation force from the waves (3.13), the drag force due to the relative velocity between the buoy and the water (3.6), the hydrostatic force (3.5) and the radiation force (3.1) from radiating waves from the interaction between the buoy and the incoming waves.

(29)

Chapter 4

Model Predictive Control of the CorPower Ocean Wave Energy Converter

Inchapter 3, the general forces acting on the WEC were described. These forces will now be simplified to form the models and cost functions that are needed for the model predictive controllers. A more exact model will choose a better control signal, but the controller will also be more computational heavy. Since this is an highly nonlinear system, two different controllers will be designed, a linear MPC and a nonlinear MPC. Both controllers will have nonlinear and nonconvex cost functions and take into consideration the constraints that exist on the system. We will assume perfect prediction of the excitation force and no measurement noise. Therefore, state estimation would be redundant.

4.1 Linear Model Predictive Control of the WEC

The first controller will be based on two time invariant linear models with a non convex cost function. As described inchapter 3, the flywheels have two different modes, engaged and disengaged. For the linear model, we will separate these two modes and create one controller for the engaged mode, and one controller for the disengaged mode.

The optimization of the torque on the generator with a disconnected flywheel will be separated from the optimization of the torque on a generator with an engaged flywheel. This is because the flywheels have non smooth and nonlinear dynamics, and it is impossible to create one single LTI system with these dynamics. The engaged MPC will predict the time until the next time each flywheel will engage, and provide this as an input for the disengaged controllers, as showed in Figure 4.1.

(30)

Figure 4.1: The connection between the engaged LMPC and the disengaged LMPCs.

4.1.1 Linear Model Predictive Control for the Engaged Mode

Model of the Engaged Mode

For the engaged mode, we will assume that there is always one generator and flywheel engaged to the rack. This engaged mode will be modelled as an linear damper acting directly on the rack of the buoy. The transmission force will be modelled as a simple linear damper and can be expressed as

Fˆ_trans= k_genv_eng, v_eng(t) ∈ V (z, ˙z, t), (4.1) where veng is the torque applied to the engaged generator. It must have the same phase as the buoy, which means that sign(veng) = sign( ˙z). If no flywheel is engaged, then veng= 0.

Other nonlinear forces acting on the buoy are the hydrostatic force, the pretension spring and the WaveSpring (F_net), as well as the hydrodynamic drag force and friction forces. We will estimate F_net as a linear equation, as shown in Figure 4.2. Note that the linear estimation ˆF_net of F_net is symmetric as it does not take the gravity into consideration. The friction will be modelled as a 5% loss of the transmission force. The drag force (3.6) will be neglected in due to it being small for low velocities.

The dynamic linear model for the engaged mode, based on (3.14), with the state variables x = (z, ˙z, q)^T, the engaged torque veng and the excitation force wexc

can now be expressed as

˙

x = Ax + B_excw_exc+ B_engv_eng (4.2)

=





0 1 0

knet

M 0 _M¹ · Cr

0 Br Ar



x +



 0

1 M

0



w_exc+



 0

kgen·kf

M

0



v_eng (4.3)

y = Cx = 0 1 0 x (4.4) subj. to: x(0) = x0 x(t) ∈ X(t), veng(t) ∈ V (x(t), t) (4.5)

(31)

Figure 4.2: Fnet with two different estimations of ˆFnet, ylinear and ypoly.

The set V (x(t), t) of possible values for v, is derived fromFigure 3.2and is V (z(t), ˙z(t), t) = {v ∈ R : |v| < vmax, sign(v) = sign( ˙z), if z ˙z ≥ 0 : v = 0} .

(4.6) We assume that there is no measurement noise in the system, and that we have a perfect prediction of the excitation force.

The continuous system is discretized using zero order hold for both the excitation force and for the torque, resulting in

x(kh + h) = Adx(kh) + Bdexcwexc(kh) + Bdengveng(kh) (4.7) where

Ad= e^Ah, Bdexc=

h

R

0

e^AsBexcds, Bdeng=

h

Z

0

e^AsBengds (4.8)

and

y(kh) = Cx(kh). (4.9)

Objective Function for the Linear Engaged MPC

For the engaged mode, the torque of the generators act as direct dampers on the rack. The overall objective of the MPC is to create a high average power output with as high quality as possible. This is done by maximizing the energy produced by the system. The total energy generated in the engaged mode can be calculated as

W =

T

Z

kgenveng(t) ˙z(t) · η( ˙ξ, veng)dt. (4.10)

(32)

Figure 4.3: Estimation of the efficiency. The point close to zero have been added to ensure that the function passes close to origo.

The efficiency η( ˙z, veng) is calculated by fitting a function of power to the graph in Figure 3.3, using the basic fitting tool in MATLAB. If we neglect low values of the rotational speed of the generator, the efficiency of the generator can be assumed to only depend on the applied torque. Figure 4.3 shows the measured efficiency values and the fitted function

ˆ

η(v(t)) = av^b(t) + c. (4.11)

Together with zero order hold for the torque v_eng(kh) of the engaged generator, and transforming the maximization problem to a minimization problem we get the discrete control problem, that optimized the average generated power,

min

veng(k) −_N¹

N

P

k=1

k_genv_eng(kh)y(kh)ˆη(v_eng(kh)) (4.12)

subj. to











x(h(k + 1)) = Adx(kh) + Bdexcw(kh) + Bdengveng(kh) y(kh) = Cx(kh)

x(0) = x0 x(kh) ∈ X v(kh) ∈ V (x(kh), t)

. (4.13)

This objective function will aim to maximize the mean produced electrical energy.

4.1.2 Linear Model Predictive Control for the Disengaged Mode

Model of the Disengaged Mode

The flywheel disconnect from the rack of the WEC when switching to the disengaged mode. In the disengaged mode the energy generation is only

Model Predictive Control of CorPower Ocean Wave Energy Converter

Model Predictive Control of CorPower Ocean Wave Energy Converter

CAROLINA ERIKSSON

Model Predictive Control of CorPower Ocean Wave Energy Converter (WEC)

Contents

List of Figures

List of Tables

Nomenclature

Abbreviations

Chapter 1

Introduction

1.1 Background

1.2 CorPower Ocean AB

1.3 Research Objectives

1.4 Thesis Overview

Chapter 2

Theoretical Background

2.1 Coordinate system

2.2 To Create Energy from Waves

2.3 Control Theory

Chapter 3

The CorPower Ocean Wave Energy Converter

3.1 Hydrodynamic Forces

3.2 Mechanical Forces of the Corpower Ocean WEC

3.3 Kinematic Relationships

Chapter 4

Model Predictive Control of the CorPower Ocean Wave Energy Converter

4.1 Linear Model Predictive Control of the WEC