Different Model Predictive Control Approaches for Controlling Point Absorber Wave Energy Converters Diploma Thesis
by
cand. kyb. Markus Richter
Examiner:
Prof. Dr.-Ing. O. Sawodny
Supervisor:
Prof. Dr.-Ing. M. E. Maga˜ na
Institute for Systemdynamics, University Stuttgart Prof. Dr.-Ing. O. Sawodny October 5, 2011
Abstract Ocean wave energy is a promising renewable energy source that can be converted into useful electrical energy using wave energy converters (WECs). In order for WECs to become a commercially viable alternative to established methods of energy generation, operating the WEC in an optimal fashion is a key task. Thus, model predictive control (MPC) is an encouraging control method that exploits the entire power potential of a WEC on the one hand, while respecting the constraints on motions and forces on the other hand. This work focuses on one subclass of WECs, namely, buoy type point absorbers with a linear generator as the power take-off (PTO) system and presents different MPCs whose goals are to optimize the generator force. For a linear one-body model of the point absorber it is shown that optimizing the power generation is possible without using an optimal velocity reference trajectory. Therefore, the control of a linear two-body model with a similar approach can be demonstrated, even though no adequate reference trajectory can be calculated in this case. Moreover, possibilities to deal with the occurance of infeasibility problems are introduced and implemented. Due to possible nonlinear effects, such as the mooring forces, a nonlinear model predictive controller is proposed, whose performance is compared to a linear MPC, also controlling the nonlinear system. The proposed controllers are validated and compared through simulation for regular and irregular sea states.
Declaration Herewith I affirm that I have written this thesis on my own. I did not enlist unlawful assistance of someone else. Cited sources of literature are perceptly marked and listed at the end of this thesis.
Stuttgart, September 2011
Acknowledgement I would like to ...
Contents Nomenclature
VII
1 Introduction
1
1.1
Ocean Energy Resources . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.2
Wave Energy Converters . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.2.1
Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.2.2
Examples of Wave Energy Converters . . . . . . . . . . . . . . .
4
Outline of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
1.3
2 Introduction to Linear Model Predictive Control 3 Linear Model Predictive Control of One-Body WEC
9 15
3.1
One-Body Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
3.2
Optimal Power Control . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
3.3
MPC Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
3.3.1
Optimal Trajectory Formulation . . . . . . . . . . . . . . . . . .
20
3.3.2
Direct Power Maximization Formulation . . . . . . . . . . . . .
23
3.4
Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
3.5
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
4 Linear Model Predictive Control of Two-Body WEC 4.1
31
Two-Body Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
4.1.1
Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . .
31
4.1.2
State Space Form and Discretization . . . . . . . . . . . . . . .
34
4.2
MPC Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
4.3
Infeasibility handling . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
4.4
Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
4.5
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
40
5 Nonlinear Model Predictive Control of a Nonlinear Two-Body WEC Model
49
II
CONTENTS 5.1
Introduction to Nonlinear Model Predictive Control
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
5.2
Nonlinear Two-Body Model . . . . . . . . . . . . . . . . . . . . . . . .
50
5.3
Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
5.4
Problem Formulation and Implementation . . . . . . . . . . . . . . . .
53
5.5
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
6 Conclusions and Future Work A Explantion: Direct Power Formulation is a Convex Problem
61 i
List of Figures 1.1
Global distribution of approximate yearly average wave power in kW/m crest length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2
3
Photograph of the Limpet [28] (a) and a schematic diagram of its principle of operation (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.3
Photographs of the Pelamis [43, 28] . . . . . . . . . . . . . . . . . . . .
5
1.4
Photograph of the Wave Dragon [28] (a) and a schematic diagram of its principle of operating (b) [17] . . . . . . . . . . . . . . . . . . . . . . .
6
1.5
Point absorber Archimedes Wave Swing [28] . . . . . . . . . . . . . . .
6
1.6
Schematic sketch (a) and photograph (b) of the L10 point absorber [28]
7
2.1
Principle of MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.2
Usual structure of MPC . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.3
Structure of MPC without reference trajectory . . . . . . . . . . . . . .
13
3.1
L10 Wave Energy Converter [38] . . . . . . . . . . . . . . . . . . . . . .
15
3.2
Electrical equivalent circuit
. . . . . . . . . . . . . . . . . . . . . . . .
18
3.3
SIMULINK simulation model for one-body direct power formulation . .
24
3.4
Simulation results with the optimal control law and a monochromatic sinusodial wave as input signal . . . . . . . . . . . . . . . . . . . . . . .
3.5
Phase relation between excitation force and actual velocity with the optimal control law . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6
27
Comparison between “MPC Traj” (T) and “MPC Power” (P) in case of monochromatic sinusodial wave as input signal . . . . . . . . . . . . . .
3.7
26
28
Comparison between “MPC Traj” (T) and “MPC Power” (P) in case of a superposition of sine waves as input signal . . . . . . . . . . . . . . .
29
4.1
L10 Wave Energy Coverter [38] (left) and schematic diagram (right) . .
32
4.2
Mooring configuration with three cables (left) and schematic diagram for force derivation (right) . . . . . . . . . . . . . . . . . . . . . . . . .
33
4.3
SIMULINK simulation model for two-body WEC control with MPC . .
39
4.4
Wave elevation, buoy position and spar position for the relaxed constraints formulation (Case 1) . . . . . . . . . . . . . . . . . . . . . . . .
42
IV
LIST OF FIGURES 4.5
Frequency of the infeasibility problem for Case 1 . . . . . . . . . . . . .
4.6
Comparison between relaxed (Case 1) and soft (Case 2) constraints formulation of the MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.7
43 44
Comparison between the soft constraints formulation with the changed parameter Thor = 6s (Case 3) and the same formulation with changed parameter ∆t = 0.2s (Case 4) . . . . . . . . . . . . . . . . . . . . . . .
4.8
44
Comparison between the soft constraints formulation with mooring constant Km = 300,000 N/m (Case 2) and the same formulation with Km = 50,000 N/m (Case 5) . . . . . . . . . . . . . . . . . . . . . . . .
4.9
45
Comparison of Spar position between the soft constraints formulation with mooring constant Km = 300,000 N/m (Case 2) and the same formulation with Km = 50,000 N/m (Case 5) . . . . . . . . . . . . . . . .
46
4.10 Comparison between the original soft constraint formulation with original WEC parameters of Table 4.1 (Case 2) and the same formulation with the changed WEC parameters of Table 4.4 for the L10 (Case 6) .
47
5.1
Comparison between the extended and the reduced model . . . . . . .
52
5.2
Nonlinear and linear mooring forces for different choices of Km . . . . .
56
5.3
Comparison between the simulation results with the NMPC and the linear MPC with Km = 100,000 N/m. The wave data is from the NDBC Umpqua buoy 46229 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4
Comparison between the simulation results with the NMPC and the linear MPC with Km = 50,000 N/m and Km = 350,000 N/m . . . . . .
5.5
57 58
Comparison between the simulation results with the NMPC and the linear MPC with Km = 150,000 N/m. The wave data is from the NDBC Umpqua buoy 46050 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6
59
Generated Power depending on Km compared to the generated power with the NMPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
List of Tables 3.1
Hydrodynamic parameter of the point absorber . . . . . . . . . . . . .
25
3.2
Simulation parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
4.1
System parameter of the two-body L10 WEC . . . . . . . . . . . . . .
40
4.2
MPC parameter values for the formulation with relaxed constraints (Case 1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41
4.3
MPC parameter values for the formulation with soft constraints (Case 2) 41
4.4
Comparison of real plant and model parameters . . . . . . . . . . . . .
47
4.5
Overview Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
5.1
NMPC/MPC parameter values . . . . . . . . . . . . . . . . . . . . . .
55
VI
LIST OF TABLES
Nomenclature
Abbreviations AWS
Archimedes Wave Swing
LWT
Linear Wave Theory
MPC
Model Predictive Control
NMPC
Nonlinear Model Predictive Control
OWC
Oscillation Water Column
PTO
Power Take-Off
R&D
Research and Development
Latin Letters Vˆ Velocity phasor F
Augmented f -vector
U
Augmented u-vector
W
Augmented w-vector
T
Augmented vector of optimal velocity
V
Augmented v-vector
X
Augmented state vector
A
Standard system state space matrix
B1
System state space matrix of input v
B2
System state space matrix of input w
B
System state space matrix of input u
A
Vector field for nonlinear discretization method
f
Input vector field nonlinear system
l
Nonlinear vector field for mooring term
p
Vector of optimization variables for NMPC
x
State vector
y
Output vector
A
Added mass (kg)
B
Damping of one-body WEC (N/(m/s))
b
Damping (N/(m/s))
VIII
Nomenclature
c
Constraint
f
System input excitation force one-body
Fe
Excitation force (N)
Fr
Radiation force (N)
fr
Radiation force impulse response function
Fgen
Generator force (N)
Fh
Hydrostatic force (N)
Fm
Mooring force (N)
g
Acceleration of gravity (m/s2 )
H
Wave height (m)
h
Discretization step time (s)
J
Objective function
j
Imaginary unit
K
Mooring cable stiffness of one cable (N/m)
k
Hydrostatic stiffness (N/m)
Km
Overall mooring cable stiffness (N/m)
L
Cable length (m)
M
Discretization order
m
Mass (kg)
me
Equivalent mass (kg)
N
Length of time horizon
Nr
Relaxation number
P
Active power (W)
Q
Reactive power (W)
q
Power, respectively, trajectory deviation weighting factor
r
Input weighting factor
ru
Relaxation factor generator force
rv
Relaxation factor velocity
rbuoy
Buoy radius (m)
S
Overall power (W)
s
Overall manipulable inputs
T
Simulation time (s)
Tp
Wave period Time (s)
Thor
Time horizon (s)
u
System input generator force
v
System input excitation force buoy
w
System input excitation force spar
w
Weighting factor slack
x
Vertical position spar (m)
Nomenclature z
Vertical position buoy (m)
Greek Letters α
Cable angle (◦ )
∆t
Optimization step time
ǫ
Slack variable
η
Wave surface elevation (m)
ω
Wave frequency (1/s)
φ
Phase angle (rad)
ρ
Density of sea water (kg/m3 )
Subscripts and Superscripts (i)
Denotes the i-th component of a vector
∗
Conjugated complex
1
Buoy parameter
12
Influence from spar on buoy
2
Spar parameter
21
Influence from buoy on spar
A
Amplitude value
d
Discrete matrices
gen
Generator
hd
hydrodynamik
k
Discrete sample steps
mcl
Meter crest length
nl
Nonlinear
opt
Optimal
p
Position
T
Transposed
u
Generator force
v
Velocity
ˆ
Denotes phasor description
IX
X
Nomenclature
Chapter 1 Introduction The idea of converting wave energy into usable electrical energy is very old. The first patent for wave energy conversion was applied in 1799 by Gilard&Son in France. But only in the 1960s, the first wave energy device, a self-recharging navigation light buoy, was built and commercialized by Masuda in Japan [31]. This device was equipped with an air turbine and today it would be considered as a (floating) oscillation water column (OWC) [14]. Since then, many patents were registered, though little attention was paid by the international scientific community for a long period of time. This changed for the first time during the oil crisis of 1973, as the interest in renewable energy technologies increased once again. After that, several research programs with governmental or private support started in Europe, Japan and the USA. One important step for wave energy conversion in Europe was taken in 1991 when the European Commission decided to include wave energy in their official research and developement (R&D) program on renewable energies [14]. Since then, the commission supported many international wave energy conferences and funded several research projects. So far, R&D in Europe is mainly being done in United Kingdom, Portugal, Ireland, Norway, Sweden and Denmark. In the USA, the Energy Act of 2005 helped support the developements in wave energy conversion, since it allowed wave energy for the first time to receive governmental support such as other renewable energy technologies [41]. Despite the large research efforts, WECs are still in the R&D stage but are closer to commercialization than ever before [12]. Thus, immense work in R&D is still required, where academic institutions make large contributions. Academic leaders in the field of wave energy conversion are currently: University of Edinburgh, Uppsala University, the Norwegian University of Science and Technology as well as Oregon State University. Since 1998, Oregon State University (OSU) in Corvallis, OR, has established a leading wave energy program for their students, where much research has been done on innovative wave energy technologies [6]. Moreover, OSU collaborates with different industrial partners such as Columbia Power Technologies. OSU also provided the
2
CHAPTER 1. INTRODUCTION
environment for creating this work.
1.1
Ocean Energy Resources
The ocean is an enormous, though so far almost untapped source of energy. Several forms of ocean energy such as tidal and marine currents, thermal, salinity and wave energy exist. This work exclusively focuses on wave energy. Wave energy is basically concentrated solar energy. Wind which creates waves by blowing across the surface of the ocean, is the result of the uneven heating of the earth’s surface. Energy conversion attempts to convert the potential and kinetic energy of the heaving water of a wave to usable energy. According to [6], the power for a monochromatic wave can be calculated by ρg 2 H 2 Tp [P/m], (1.1) Pmcl = 32π where the wave height H is measured from trough to crest, ρ is the density of sea water, g is the accaleration of gravity and Tp is the wave period. The power is usually given in the unit of Power per meter crest length and is proportional to wave height squared and linearly dependent on the wave period. It is difficult to estimate the amount of exploitable wave energy in the world’s ocean. According to [25], the ocean holds approximately 8, 000 − 80, 000 TWh/year or 1 − 10
TW, whereas Falnes in [32] quantifies the world’s exploitable wave power resource to be of the order of 1 TW. Comparing this to the world’s annual energy consumption of approximately 148, 000 TWh in 2008 [40] shows that wave energy could play an important role in the world’s energy portfolio. In comparison to solar and wind energy, wave energy stands out due to its low hourly variation, high power density and its high availablity [6, 23]. Additionally, wave’s motion can be forcasted very accurately. It should also be noted that average
wave power is cyclical, where the waves in the winter are larger than in the summer. Moreover, the wave power decreases as the waves advance toward the shoreline on account of frictional losses with the sea floor [25]. Due to the west-to-east winds, the waves on the west coasts of the continents contain more energy as can be seen in Fig.1.1. Oregon with a yearly average power of about 40 kW/m is characterized by a similar average power level as, for instance, Portugal or Norway. The British Isles, for example, provide an even higher power level of about 70 kW/m. The high power levels on these coasts may be a reason why wave energy research has especially developed in these countries.
1.2. WAVE ENERGY CONVERTERS
3
Fig. 1.1: Global distribution of approximate yearly average wave power in kW/m crest length [29].
1.2
Wave Energy Converters
In order to generate power from wind, the wind turbine is the only commercially viable design. The situation for wave energy conversion is very different. There are various concepts which differ in their power conversion principle, in their size or in their location on the ocean. However, there is no concept which has really come out on top so far. In the following, WECs are classified and currently promising WECs are described.
1.2.1
Classification
Converters are generally classified according to their location with respect to the shore or to the applied conversion principle. Both ways of classification are discussed here following the work in [31]. Shoreline devices can be located directly on the sea floor, in shallow water or can be attached to a rock. Thus, they are easy to install and to maintain. Additionally, short cables can be used and the environment for the device is very safe. As stated above, waves loose energy when advancing toward the shoreline. Hence, the environment for generating power near the shore is generally weak. However, in some areas the coastal geometry helps focus the energy and thereby compensates for the weak environment. Other disadvantages lie in the difficulty of finding suitable sites and in the impact on coastal landscapes and environments. Offshore devices are floating on the sea or are submerged in deep waters and need to be moored to the sea floor. They can harness the high power potential of the open
4
CHAPTER 1. INTRODUCTION
sea and do not have an effect on the coastal landscapes and environments. However, installation and maintainance costs are very high and harsh sea conditions give rise to survivability and reliability problems. According to [6], four different operating principles for WECs exist: • Oscillating Water Column (OWC)
The water level in a chamber rises and falls with the incoming waves. Air rushes in and out a small opening in the chamber driving a turbine which generates electricity.
• Attenuator
It is a device with floaters, where the motion of the floaters relative to each other causes energy generation by hydraulic pumps or other type of converters. It is usually a rectangular device which can be oriented perpendicular or colinear with the wave front.
• Overtopping
The device concentrates water with a ramp in a higher reservoir and releases the water through a hydro turbine.
• Point Absorber
A relatively small device which captures energy from all directions at one point by a floater which is moved by the waves. Usually the motion is only heave motion.
1.2.2
Examples of Wave Energy Converters
In what follows, some WECs are briefly described. The first commercial scale wave energy converter was the Limpet (Land Installed Marine Powered Energy Transformer) [42] which is located on the shoreline of the Island of Islay, Scotland. The device was built in 2000 and has produced power for the national grid since then. Its nominal power is 500 kW. Limpet belongs to the class of OWC devices. Fig. 1.2 shows a picture of the Limpet and a schematic diagram of its principle of operation. The water level in the chamber rises and falls with the incoming waves causing air compressing and decompressing which drives a turbine. Since the direction of the air flow changes permanently, a special type of turbine, namely the Well’s turbine, is used to assure only one direction of rotation. Limpet is also used as a testing facility for the development of other turbines. The Pelamis shown in Fig. 1.3 is floating offshore. It is an attenuator device with four segments joined together by hinged joints [30]. By motion of the waves, the joints flex relative to each other. The motion in the joints is resisted by hydraulic cylinders
5
1.2. WAVE ENERGY CONVERTERS
(a)
(b)
Fig. 1.2: Photograph of the Limpet [28] (a) and a schematic diagram of its principle of operation (b).
(a)
(b)
Fig. 1.3: Photographs of the Pelamis [43, 28].
that pump fluid into high pressure accumulators. By means of a hydraulic engine power can be generated [30]. The first full-scale prototype was completed 2004 and deployed in Orkney, Scotland with the rated power of 750 kW. The world’s first small wave farm with three devices was deployed in Portugal 2008. Additionally, further projects, such as the wave farm with 14 devices in Shetland, Scotland, are in the development stage. The Wave Dragon, Fig. 1.4, represents an overtopping device. The first prototype was deployed at Nissum Bredning, Denmark in 2003. The device is floating on the sea. By means of two “arms” (reflectors) and a submerged ramp, water is guided up into a reservoir. The reservoir is higher than the sea level, thus potential energy is stored which can be transformed into usable energy by letting the water run through the turbine outlet. Further projects are planned in Wales and Portugal. The Archimedes Wave Swing (AWS) belongs to the class of point absorbers. Two photographs can be seen in Fig. 1.5. It is a submerged device and was first deployed in Portugal in 2004 and successfully connected to the grid shortly afterwards. The
6
CHAPTER 1. INTRODUCTION
(a)
(b)
Fig. 1.4: Photograph of the Wave Dragon [28] (a) and a schematic diagram of its principle of operating (b) [17]
(a) Several AWS under water.
(b) AWS plant before submersion.
Fig. 1.5: Point absorber Archimedes Wave Swing [28].
operating principle is based on pressure differences in a chamber when the wave rises and falls. The chamber of the AWS is filled with air, where the bottom is fixed and the lid moves up and down due to the pressure difference on top [31]. In the first prototype the PTO system was a linear generator. Later, hydraulics were used as well. One advantage of this type of point absorber is its completely submerged design which is in general less vulnerable in harsh sea conditions and does not have an negative impact on the landscape. Fig. 1.6 shows the point absorber L10, developed at Oregon State University in collaboration with Columbia Power Technologies and the Navy, where a linear generator as the direct-drive PTO system is used. This device was tested in September 2008 off Newport, OR. Since then, Oregon State University is continuing the collaboration in order to develope a full-scale WEC [6]. Various WECs have been deployed and tested in the last 10 years. However, there is no concept or design which came out on top so far. Moreover, the design of wave
7
1.3. OUTLINE OF THE THESIS
(a)
(b)
Fig. 1.6: Schematic sketch (a) and photograph (b) of the L10 point absorber [28]
farms is still in its early stages. It will be interesting to observe which concept will be the most successful one in the future. Even though this work exclusively focuses on the L10 point absorber, certain principles could also be applied to other designs.
1.3
Outline of the Thesis
Basically, controllers for WECs need to perform two main tasks, i.e., the generated power needs to be maximized and the physical constraints on the machinery force and on the WEC’s motion need to be respected. Additionally, controllers which can benefit from prediction data of the wave motion are promising, since wave prediction is generally possible. For this reason, MPC is a reasonable type of control algorithm, since it provides a framework for optimal control and constraints. Moreover, prediction data can be naturally included. The objective of this work is to implement MPCs for controlling a point absorber which is capable of operating it under optimal conditions, where the power generation is maximized, while the constraints of the WEC are satisfied. Depending on the model of the device, the MPCs differ. The main part of this thesis is devided into four parts. The basic principles and strategies of model predictive control (MPC) are briefly discussed in the beginning. Chapter 3 discusses exclusively a linear one-body model of a point absorber WEC. Based on the one-body model, an optimal control law for power optimization is introduced and two different MPC formulations are implemented. The focus lies on demonstrating that MPC of the WEC is possible without using an optimal velocity trajectory as reference. The different controllers are validated and compared by means of simulation results with a superposition of sine waves. Based on linear wave theory (LWT), a linear two-body model for a point absorber
8
CHAPTER 1. INTRODUCTION
with sea floor mooring is derived and a MPC is implemented in Chapter 4. In MPC applications, infeasibility problems can occur while solving the optimization problem. Since this can cause big problems, two possibilities to handle the occurance of infeasibilities are proposed and compared by simulation results for irregular sea. Additionally, the influence of MPC parameters, mooring and model inaccuracies is discussed. In Chapter 5, due to possible unmodelled nonlinear effects, a nonlinear model is derived, where the mooring and radiation forces are assumed to be nonlinear. In order to control the nonlinear system, two different controllers are proposed and compared. First, the linear MPC from Chapter 4 and secondly a nonlinear MPC whose implementation is the focus in this chapter. By variation of the mooring constant of the linear MPC, its performance is analyzed and compared to the nonlinear model predictive controller (NMPC). Chapter 6 concludes this thesis and suggests areas for future research.
Chapter 2 Introduction to Linear Model Predictive Control The first theories of model predictive control were introduced in the 70’s and were also known as receding horizon control and open loop feedback optimization [3]. At that time, most of the research was done in the area of process control engineering. Since then, the theories have advanced and the computational power of control units has increased tremendously so that today’s MPC applications are used in various technical control problems. One main reason for the extensive capabilities of MPC is its flexibility, namely, • System states constraints or actuator limits can be taken into account • It can handle multivariable control • It can handle input and output time delays • It can handle unstable and non-minimal phase processes • It can include predictions of external inputs naturally Despite its flexibility, the core compenents of MPC are always the same, i.e., • Prediction model of the process • Objective function • Optimizing algorithm In fact, many different MPC algorithms exist, where the model representation and the form of the objective function, as well as the optimizing algorithm can vary. Nevertheless, all MPC algorithms have in common the same control idea. The current control signal is obtained by solving a finite horizon open loop optimal control
10
CHAPTER 2. LINEAR MPC
problem (mostly, minimizing an objective function) online at each time step [9]. Hence, the current system state is used as the initial state for the model. This yields a sequence of optimal control signals at every time step, where only the first signal is applied to the process. Also, the prediction horizon is shifted one time step before the new optimization starts. Fig. 2.1 shows this basic principle of MPC. In most cases, the u(t+k∆T |t)
u(t)
w(t+k∆T ) y(t) y(t+k∆T |t)
t−∆T
t
t+∆T . . .
N
...
t+k∆T
t+N ∆T
Fig. 2.1: Principle of MPC.
goal is that the system outputs y(t + k∆T | t)1 for k = 1, . . . , N follow the predefined
reference trajectory w(t + k∆T ) for k = 1, . . . , N . N denotes the number of time steps within the prediction horizon and ∆t the optimization step time. A typical objective function consists of two quadratic terms. One term evaluating the reference deviation and the other term weighting the energy use. A possible objective function for a single-input single-output system, similar to the one in [9], is given by J(N ) =
N X j=1
N 2 X 2 δ(j) y(t + j∆T | t) − w(t + j∆T ) + λ(j) ∆u(t + (j − 1)∆T ) , (2.1) j=1
where the coefficients δ(j) and λ(j) are sequences of mostly constant numbers weighting the reference deviation and the energy use, respectively. If the objective function is quadratic and there are no constraints, the control signals can be calculated explicitely [9]. Otherwise, iterative optimization algorithms have to be used. Other types of objectice functions can be found in [3, 8], where an explicit terminal point condition or separate control horizons for prediction and controlling are included. In order to handle, for instance state transitions, the reference trajectory can be just a step or a smooth approximation of that. Other differences between MPC algorithms exist in the description of the model. In general, most descriptions can be used, but it is of great importance to have an appropriate model description in order to describe the dynamic behaviour of the pro1
This notation denotes the variable values at the time t + k∆T , calculated at the time t
11 cess adequately. In order to model the input/output behaviour of a process, system identification by means of an impulse or a step response is often used. Also, transfer functions are often used. As is mostly the case, state space models are applied to characterize the system. A possible form for a continous, linear, time-invariant system is x˙ = Ax + Bu,
x(0) = x0 , x ∈ Rn ,
y = Cx,
(2.2) u ∈ Rm ,
y ∈ Rr .
(2.3)
In practice, there are constraints for most processes as well. Actuators have certain physical or safety limits and also the states are often restricted, so that in most cases constraints need to be considered in the optimization algorithm. Normally, the constraints are state and input constraints as umin ≤
∀t,
(2.4)
∆umin ≤ ∆u(t) ≤ ∆umax ∀t,
(2.5)
xmin ≤
u(t) ≤ umax x(t) ≤ xmax
∀t.
(2.6)
The introduction of constraints yields more difficult optimization problems which are only solvable with iterative methods, such as the active-set or the interior-point method, as described in [27]. Now that the principles of MPC have been discussed, it has to be stated why MPC is an adequate and promising control algorithm, especially in the field of point absorber WECs, [7]. • The main goal of controlling point absorber is to maximize the power generation
in order to make WECs cost-effictive and able to compete with other methods of energy conversion. At the same time, mechanical limits for PTO force and system states need to be considered. For that, MPC provides a framework to explicitely consider these constraints while maximizing the power generation.
• In general, having a prediction of the incoming waves yields benefits for optimal energy capturing. However, many control algorithms are not capable of exploiting
prediction information. As stated above, predictions of external inputs can be easily included in the MPC framework though. • Traditionally, MPC has been applied to chemical or process control engineering applications, since those are slower dynamic processes which do not overburden
the computational power for online optimization. Even though faster processes are controlled with MPC nowadays, it is always easier to apply MPC if the process is not too fast. Ocean wave applications are in general rather slow processes and therefore fit well the MPC requirements.
12
CHAPTER 2. LINEAR MPC
In the following chapter, two different versions of MPC are introduced. As an outlook, the two different MPC structures are briefly discussed here. The first approach requires a reference, namely an optimal velocity reference trajectory which the actual velocity has to follow in order to achieve optimal energy capturing. The calculation of this optimal trajectory is discussed in Section 3.2. In this case, the usual and above described MPC structure can be used which is schematically shown in Fig. 2.2. The objective Current States
Predicted Disturbances
Reference Trajectory
Model
Future Outputs
Future Inputs Optimizer
+
−
Error
Optimal Input
Objective Function
Constraints
Fig. 2.2: Usual structure of MPC.
function contains a term penalizing the error between reference trajectory and actual velocity. The optimizer tries to find the optimum of the objective function with respect to the constraints. Also, the MPC gets the current state information from measurements and calculates with the aid of the process model the future outputs which are compared to the reference. After optimization, the first control signal is applied to the process. This work shows that even in the case when no optimal trajectory can be calculated, MPC can still be applied. Therefore, the objective function needs to be changed so that one term directly indicates the generated power, as can be seen in Section 3.3.2 and thereafter. The structure of this second version is shown in Fig. 2.3, where no reference is used so that only the input and the system states are part of the objective function. Here, future information is used as well. However, not by a reference trajectory but only by the prediction of external disturbances which is also discussed in Section 3.3.2.
13
Current States Predicted Disturbances
Model
Future States
Future Inputs Optimizer Optimal Input
Objective Function
Constraints
Fig. 2.3: Structure of MPC without reference trajectory.
14
CHAPTER 2. LINEAR MPC
Chapter 3 Linear Model Predictive Control of One-Body WEC This chapter mainly discusses the application of linear MPC to a relatively simple onebody representation of a point absorber. At first, a state space model is derived and brought into a discrete representation before different control methods are proposed. The following two MPC formulations specifically differ in the use of prediction data and in the formulation of the objective function. Also, the MPCs are compared with each other and to an optimal control law via simulation.
3.1
One-Body Model
Fig. 3.1 shows a drawing of the point absorber L10, developed at Oregon State University. The L10 consists of two bodies, a buoy floating on the surface of the ocean
Fig. 3.1: L10 Wave Energy Converter [38]
16
CHAPTER 3. ONE-BODY WEC
and a damping body, namely, a combination of a spar and a ballast tank. The two bodies are connected through a PTO system in order to convert the relative motion into usable energy. The L10 is an example of a direct-drive WEC. In this case, the wave energy is converted directly into electrical energy by means of a linear generator without converting wave energy into mechanical energy as in hydraulic PTO systems first. In what follows, the point absorber is simply considered as an axis-symmetric float connected to a stationary reference body through an electric linear generator, as considered in [7]. This is also comparable to the assumption of perfect mooring of the spar as the reference body. Additionally, linear hydrodynamics and frequency independent WEC parameters are assumed for modelling. These simplifications are reasonable for simple body shapes and small ranges of motion. Based on Newton’s law, the dynamic equations can be written as m¨ z (t) = Fe + Fgen + Fr + Fh ,
(3.1)
where z¨(t) is the buoy acceleration, m is the mass of the device and Fgen is the force produced by the power take-off system. Fgen presents also the manipulable input to control the system. According to [15], at least three other forces act on the body: The radiation force, created by the moving of the float and thus radiating waves; the hydrodynamic force which represents the restoring force of the water; and the excitation force, which is the force the incoming wave exerts on the body. The excitation force Fe represents an unmanipulable system disturbance. The radiation force Fr and the hydrodynamic force Fh can be written simply as Fr = −A¨ z (t) + B z(t), ˙
(3.2)
Fh = −gρπrf2 loat z(t) = −kz(t),
(3.3)
where A is the added mass of the body, B is the viscous damping, g is the acceleration of gravity, ρ is the density of seawater, rf loat is the radius of the float and k is then called the hydrostatic stiffness. By means of the Morison approach [15], the wave motion can be linearized and the excitation force can be described by Fe = A¨ η (t) + B η(t) ˙ + kη(t),
(3.4)
where η is the wave elevation. Summarized can be obtained Fe + Fgen = (A + m)¨ z + B z˙ + kz.
(3.5)
The system described by (3.5) can be written in state-space form as " # " #" # " # # " z˙ 0 1 z 0 0 = + Fgen + Fe . k B 1 1 z¨ − m+A − m+A z˙ | | m+A | m+A {z } {z } {z } A
B
B
(3.6)
17
3.2. OPTIMAL POWER CONTROL
The output of the system is equal to the states of the system " # z y=x= . (3.7) z˙ The MPC needs a discrete system representation whose sampling time h has to be the same as the optimization sample time ∆t for the MPC. In what follows, the inputs Fgen and Fe at the sample instant k are denoted by u and f , respectively. Let xk = x(kh). Assuming a Zero-Order Hold for all inputs u(t) = uk , v(t) = vk ,
(3.8)
with kh ≤ t ≤ (k + 1)h, the following discrete system can be obtained [11] xk+1 = Ad xk + Bd uk + Bd fk ,
x0 ∈ R2 ,
(3.9)
where Ad = eAh , Z h eAφ dφ B, Bd =
(3.10)
0
Many different ways of computing the matrix exponential eAh can be found in [24]. It should also be noted that the discretization is exact for the input u, since it will be calculated with the MPC at every sample time and is therefore a piecewise constant signal as assumed in (3.8). For the excitation force f , which is in general not piecewise constant, a small error can be expected which grows with larger sample times.
3.2
Optimal Power Control
In the past much work has been done on control algorithms, such as latching control for WECs, and on optimizing the energy absorption [4, 18, 15]. Most of these studies assume, similar to here, one wave-absorbing semisubmerged heaving body. Often mooring forces are not taken into account as well. For these assumptions, an optimal control law can be derived. The next explanations follow the work in [6, 15]. Making the assumption that the wave climate is monochromatic and there is a steady state wave input, phasors can be used, since all equations are linear. Then, the position and the excitation force can be written as ηA jφη ˆ =√ (3.11) e , ⇒ H 2 Fe,A Fe (t) = Fe,A cos(ωt + φFe ) ⇒ Fˆe = √ ejφFe , (3.12) 2 , and ω and Tp denote the wave frequency and time period, respectively. where ω = 2π Tp η(t) = ηA cos(ωt + φη )
Also, φ is the phase angle and the subscript A denotes peak amplitude values. Thus,
18
CHAPTER 3. ONE-BODY WEC
(3.5) can be written in phasor form ˆ Fˆe + Fˆgen = (jω(A + m) + B + k/(jw))jω H = (jω(A + m) + B + k/(jw)) Vˆ , {z } |
(3.13) (3.14)
Zhd
where Vˆ is the phasor of the velocity η˙ and Zhd can be considered as hydrodynamic impedance, analogous to an electric circuit. In order to find an optimal control law, the linear generator is assumed to be equivalent to a spring-mass-damper system. Hence, the generator produces a force proportional to acceleration, velocity and position of the buoy. Fˆgen = − (jωmgen + Bgen + kgen /(jw)) Vˆ . {z } |
(3.15)
Zhd,gen
The minus sign is added in order to consider the generator as a resistance and not as as source. Thus, (3.14) and (3.15) can be interpreted as the equivalent of the electric circuit shown in Fig. 3.2. Voltage is equivalent to force and current is equivalent to the
Fig. 3.2: Electrical equivalent circuit
buoy velocity. Different types of power components need to be distinguished. S = Fˆe Vˆ ∗
(3.16)
Shd = Phd + jQhd = (jω(A + m) + B + k/(jw))Vˆ · Vˆ ∗ = (B + j(ω(A + m) − k/ω))|Vˆ |2 ,
(3.17)
= (Bgen + j(ωmgen − kgen /ω))|Vˆ |2 ,
(3.18)
Sgen = Pgen + jQgen = (jωmgen + Bgen + kgen /(jw))Vˆ · Vˆ ∗ where S is overall power delivered to the WEC, Shd is the power delivered to the water and Sgen is the power fed into the electrical generator. P and Q denote the active and reactive power, respectively. If the imaginary component of the generator power
3.2. OPTIMAL POWER CONTROL
19
is nonzero, the generator is delivering power to the sea at some points. Thus, external energy needs to be fed into the generator in order to source the reactive components for several seconds. Therefore, the storage of energy is necessary which can be a significant requirement for the WEC. Maximizing power generation for the equivalent circuit is an impedance matching problem, [2]. With Fˆe Fˆe , = Vˆ = Zhd + Zhd,gen (jω(A + m) + B + k/(jw)) + (jωmgen + Bgen + kgen /(jw)) (3.19) it follows that the active generator power is given by 2 ˆ Fe 2 ˆ Pgen = Bgen |V | = Bgen . Bgen + j(ωmgen − kgen /ω) + B + j(ω(A + m) − k/ω) (3.20)
By the following appropiate choice, the reactive components can be cancelled and the power generation is maximized. mgen = −(A + m),
(3.21)
Bgen = B,
(3.22)
kgen = −k.
(3.23)
For the generated power yields Fˆ 2 |Fˆ |2 e e (3.24) Pgen = B = 2B 4B and hence for the optimal velocity Fˆe Vˆ = . (3.25) 2B This means that the control law alters the resonant frequency of the system to resonate with the dominant wave frequency, ω. This control law is also called reactive control or phase and amplitude control [15], because condition (3.25) means that the velocity Vˆ (or z˙com ) has to be in phase with the excitation force and the velocity amplitude |Vˆ |
has to equal
Fˆe , 2B
[15].
Finally, the optimal control law can be written as Fgen,opt = (A + m)¨ z − B z˙ + kz.
(3.26)
It can be observed that the theoretical optimal generated power is exactly one half of the incoming power [6], where the other half diffuses back into the sea. In the case of allowing surge motion, the theoretical limit of power generation is 1. Also, it has to be mentioned that this control law requires significant energy storage because of the reactive components.
20
CHAPTER 3. ONE-BODY WEC
3.3
MPC Formulation
In this section, two different MPC approaches are proposed. The first formulation requires an optimal reference trajectory for the buoy velocity. This reference can be calculated by using (3.25). Furthermore, constraints for the generator force, buoy position and buoy velocity are considered. The second formulation attempts to optimize the power generation by directly including a power term into the objective function, and thus no reference is neccessary.
3.3.1
Optimal Trajectory Formulation
The objective function in this formulation needs to contain one term which evaluates the deviation of the actual buoy velocity from the optimal reference velocity and another term which evaluates the energy use. Additionally, constraints need to be considered. In the following, Fgen and Fe are denoted by u and f , respectively. Hence, the optimization problen can be defined as min J(xk , uk ),
xk ,uk
(3.27)
where N
J(xk , uk ) = subject to
1X (xk − xk,opt )T Q(xk − xk,opt ) + r u2k−1 , 2 k=1
(3.28)
xk+1 = Ad xk + Bd uk + Bd fk
(3.29)
|zk | ≤ cp , ∀k = 1, . . . , N
(3.30)
|z˙k | ≤ cv , ∀k = 1, . . . , N
(3.31)
|uk | ≤ cu , ∀k = 0, . . . , N − 1,
(3.32)
where Q=
"
0 0
#
, (3.33) 0 q with weighting factors q ≥ 0, r > 0 and the symmetric relative position cp , velocity cv
and input cu constraint. N denotes the length of the time horizon. xk,opt denotes the
optimal state/output trajectory which the actual state/output has to follow. Only the velocity values of the reference trajectory are weighted, since there is no information about the optimal position trajectory. So, all optimal position values can be set to zero. To reformulate the optimization problem statement, define the following four aug-
21
3.3. MPC FORMULATION mented vectors: X =
h
xT1 , xT2 , ..., xTN
h
iT
, iT
U = u0 , u1 , ..., uN −1 , h iT F = f0 , f1 , ..., fN −1 , iT h T T T T = t1 , t2 , ..., tN ,
(3.34)
with dim(X , T ) = 2N and dim(U, F) = N , where T is the optimal output state
trajectory over the entire horizon with
tk =
"
0 1 f 2B k
#
.
(3.35)
Thus, the objective function (3.28) can be formulated in matrix form as 1 ˆ (X − T ) + 1 U T RU, ˆ J = (X − T )T Q 2 2 ˆ und R ˆ are diagonal matrices where Q ˆ = diag(Q, Q, . . . , Q), Q
ˆ = diag(r, r, . . . , r) R
(3.36)
(3.37)
ˆ = (2N × 2N ) and dim(R)(N ˆ ˆ has to be at least positive semidefwith dim(Q) × N ). Q ˆ positive definite. inite and R In order to decrease the computational effort, the problem will be formulated only in terms of the input variables, similar to the work in [44, 35]. A closed-form expression for all states xk that only depends on uk and x0 can be formulated by solving the system and substituting in the system equation (3.29) x 1 = A d x 0 + B d u0 + B d f 0 , x 2 = A d x 1 + B d u1 + B d f 1 = A d 2 x 0 + A d B d u 0 + B d u1 + A d B d f 0 + B d f 1 , .. .
(3.38)
xN = Ad N x0 + Ad N −1 Bd u0 + Ad N −2 Bd u1 + Ad N −1 Bd f0 + Ad N −2 Bd f1 + . . . + Bd uN −1 + Bd fN −1 . In vector form,
x1
Ad
Bd 0 0 2 x 2 Ad Ad B d Bd 0 x A 3 A 2B Ad B d Bd 3 = d x0 + d d . . . .. .. .. xN Ad N Ad N −1 Bd Ad N −2 Bd {z | {z } | {z } | X Jx Ju
... ... ... ... ...
f 0 u 0 0 u f 0 1 1 0 u2 − f 2 . .. .. . . . . fN −1 Bd uN −1 } | {z } | {z } U F (3.39)
22
CHAPTER 3. ONE-BODY WEC
with dim(x0 ) = 2,
dim(X ) = 2N,
dim(Jx ) = (2N × 2),
dim(Ju ) = (2N × N ).
dim(U) = N,
(3.40) (3.41)
By plugging (3.39) in the objective function (3.36), X is eliminated and one gets 1 ˆ ˆ x x0 + Ju U + Ju F − T ) + 1 U T RU. (3.42) J = (Jx x0 + Ju U + Ju F − T )T Q(J 2 2 By expanding and neglecting the terms without dependence on U, the objective func-
tion can be written as 1 ˆ u + R)U ˆ + U T (Ju T QJ ˆ x x0 − Ju T QT ˆ + Ju T QJ ˆ u F), J = U T (Ju T QJ 2 and is now only quadratic dependent on U.
(3.43)
The constraints need to be reformulated as well. The constraints from (3.30) -
(3.32) can be also written as Dk uk ≤ dk , for k = 0, . . . , N − 1,
(3.44)
Ek xk ≤ ek , for k = 1, . . . , N.
(3.45)
The input constraints are already only dependent on U and can be written in matrix
form as
DD U ≤ d,
(3.46)
where DD = diag(D0 , . . . , DN−1 ) and d = [dT0 , . . . , dTN −1 ]T . The state constraints contain the state xk . To get rid of it, xk is replaced by (3.39). It follows E 1 B d u0 ≤ e 1 − E 1 A d x 0 − E 1 B d f 0 , .. . EN (Ad N −1 Bd u0 + Ad N −2 Bd u1 + . . . + Bd uN −1 ) ≤ eN − EN Ad N x0 − EN (Ad N −1 Bd f0 + Ad N −2 Bd f1 + . . . + Bd fN −1 ).
(3.47) With the definitions ED = diag(E1 , . . . , EN ) and e = [eT1 , . . . , eTN ]T the constraints can be written as ED Ju U ≤ e − ED Jx x0 − ED Ju F.
(3.48)
Now, the MPC problem can be restated as min J(U), U
(3.49)
where 1 ˆ + Ju T QJ ˆ u F), ˆ u + R)U ˆ + U T (Ju T QJ ˆ x x0 − Ju T QT J = U T (Ju T QJ 2 subject to " # # " DD d . U≤ E D Ju e − E D Jx x 0 − E D Ju F
(3.50)
(3.51)
23
3.3. MPC FORMULATION
Problems in this form are usually called Linear Quadratic Problems (LQP), because the objective function is quadratic and all constraints are linear. There exist many appropriate methods to solve problems of this class. The reader is referred to [27] and [34].
3.3.2
Direct Power Maximization Formulation
The next MPC problem formulation is also done in terms of the input variables U,
which allows the constraints formulation to remain the same. The difference is that this formulation can maximize the power directly, without the need of an optimimal trajectory for the buoy velocity. This can be very advantageous in case no optimal trajectory can be calculated. Also, the assumption of a sinusoidal monochromatic wave in order to calculate the reference can be a problem in the previous formulation. For irregular waves, it generally cannot be expected that the calculated optimal trajectory be the optimal one. Moreover, this trajectory can be far away from the optimal power generation under certain conditions. The following formulation focuses on optimizing the generated power without using a reference. In general it is defined by Pgen = −Fgen (t)z(t). ˙
(3.52)
In matrix form over the complete horizon, the power is given by T Pgen = −X(2) U
where the subscript
(2)
(3.53)
denotes the second component of each xk in X . To extract the
second components of X , the matrix S is defined, i.e. 0 1 0 0 ... 0 0 0 0 1 . . . 0 S= .. .. , . .
(3.54)
0 0 0 0 ... 1
with dim(S) = (N × 2N ), so that the following holds: Pgen = −(SX )T U.
(3.55)
Plugging (3.55) in (3.39) yields Pgen = −(S(Jx x0 + Ju U + Ju F))T U
(3.56) = −x0 T Jx T ST U − U T Ju T ST U − F T Ju T ST U = − U T Ju T ST U + x0 T Jx T ST + F T Ju T ST U . To make this problem fit the minimization problem, the sign is switched. Hence the MPC formulation for generating power optimization is min J(U), U
(3.57)
24
CHAPTER 3. ONE-BODY WEC
where 1 ˆ U + 1 x0 T Jx T ST + F T Ju T ST U, J(U) = U T Ju T ST + R 2 2
(3.58) (3.59)
subject to "
DD E D Ju
#
U≤
"
d
e − E D Jx x 0 − E D Ju F
#
,
(3.60)
ˆ is again the weighting matrix for the input. where R
The trajectory formulation is naturally convex, thus every local solution is the global solution of the optimization problem. That is why convexity is a requirement for most global optimization algorithm for LQP problems. It can also be shown that the power formulation yields a convex problem which is no longer symmetric, see Appendix A. However, this is only shown for the two-body case due to redundancy. Hence, the same optimization methods used in the trajectory formulation can be applied.
3.4
Implementation
The simulation is implemented using Matlab/SIMULINK. The simulation model for the one-body direct formulation is shown in Fig. 3.3. It mainly consists of four different (0.1
0.2
U Y
pos_wave
U Y
vel _wave
0.3)
F_gen
Amplitude
[1/(2*pi) 1/(2.5*pi) 1/(3*pi)] Frequency
(−pi/2
generate _sinus
eta
F_e
Exitflag Optimization
U Y
0 pi/2)
F_e
pos_buoy
MPC
F_e computation
Phase Iteratations first optimization
0 Offset
vel _buoy S−Function
Level −2 M−file S−Function F_gen F_e
Fgen
Fe
pos
pos_buoy
vel
vel _buoy
Exitflag Infeasibility catching
RealPlant
Fig. 3.3: SIMULINK simulation model for one-body direct power formulation.
blocks. The “generate sinus” block is able to generate a superposition of sine waves as wave input signal with corresponding prediction over a certain horizon. The wave signal consists of position, velocity and acceleration, and perfect prediction is assumed. According to (3.4), a prediction for the excitation force can be calculated in the block “Fe computation”. The actual MPC algorithm is implemented in the Level-2 M-file s-function “MPC”. As stated earlier, the inputs of the “MPC” are the predictions of the excitation force, as well as the current system state, which is the output of the “RealPlant” block. It is assumed that the current system state can be measured. The
25
3.5. RESULTS
ouputs of the “MPC” block are the first generator force signal and certain signals for the optimization status. The hydrodynamic parameters of the point absorber used in the simulation are given in Table 3.1. The parameters are taken from [7] and were determined from analysis Variable
Explanation
Value
m
Mass
1997 kg
A
Added mass
5660 kg
B
Viscous damping
11,400 N/(m/s)
k
Hydrostatic stiffness 88,970 N/m
Tab. 3.1: Hydrodynamic parameter of the point absorber.
with ANSYS AQUA [1]. The simulation parameter used for the MPC simulations are given in Table 3.2. The position constraint is equivalent to the stroke length of the Variable
Explanation
Value
∆t
Sample time for optimization
0.1 s
Thor
Optimization horizon
3s
N
Number of values
30
T
Simulation Time
50 s
cp
Position constraint
1m
cv
Velocity constraint
1 m/s
cu
Generator force constraint
80,000 N
Tab. 3.2: Simulation parameter.
generator and is taken from [33]. The generator constraint is chosen so that the MPC control can be effectively tested.
3.5
Results
In this section, the two different implementations are verified and compared by means of computer simulation. In the following, “MPC Traj” and “MPC Power” denote the MPC with the trajectory formulation and the direct power formulation, respectively. In all figures they are abbreviated by (T) and (P), respectively, and the optimal control law by (O). First, a simple sinosodial wave with an amplitude of 0.6 m and a period time of 7 s is used as the input wave signal. In Fig. 3.4, the simulation results with the optimal control can be seen. On the top left, the wave elevation η and the actual buoy displacement with optimal control are plotted. It can be noted that displacements of 2.4 m are achieved which would be clearly above the actual requirements. The
26
CHAPTER 3. ONE-BODY WEC Displacement of Wave and Buoy
5
3
3
η z (O)
1 0 −1 −2 −3 0
Excitation and Generator Force Fe Fgen (O)
2
Force [N]
Position [m]
2
x 10
1 0 −1 −2
10
20
30
40
−3 0
50
10
20
30
Time [s]
Time [s]
Wave, Buoy and Commanded Velocity
Power Generation
40
50
3
η˙ z˙ (O) z˙ com
Pgen (O)
300 200
1
Power [kW]
Velocity [m/s]
2
0 −1
100 0 −100
−2 −200 −3 0
10
20
30
40
50
0
10
Time [s]
20
30
40
50
Time [s]
Fig. 3.4: Simulation results with the optimal control law and a monochromatic sinusodial wave as input signal.
excitation force, seen on the top right, is in the range of ±50,000 N, whereas the
generator force with about 200,000 N would violate the constraints. On the bottom
left it can be seen that the actual velocity with this control law is exactly the same as the earlier calculated commanded velocity which is later used as a reference. The generated power, plotted on the bottom right, is within −210 kW and +260 kW and
the average generated power is 25.24 kW. It has to be noted that the amount of reactive power at some points is enormous, so that a large energy storage is necessary. Fig. 3.5 shows the phase relation between excitation force and actual velocity with
the optimal control law. It can be seen that the optimal velocity has to be in phase with the excitation force, which was discussed in Section 3.2. The next simulation compares the two different MPCs, where the same wave data as in the previous simulation is used. The parameters and constraints of Table 3.2 are applied, and Fig. 3.6 shows the results. The buoy positions for both controllers satisfy the limits of ± 1 m well. Also, almost no difference between both trajectories can
be seen. The same holds for the generator force. Both trajectories are almost the same and respect the generator limit as desired. On the bottom left, the buoy and the commanded velocity are plotted. The MPC attempts to follow the commanded velocity, but is restricted by the constraints. Even though, “MPC Power” is not designed to
follow the comanded velocity, the same results as for “MPC Traj” are obtained, and
27
3.5. RESULTS 2.5
Fe z˙ (O)
Force [N/105] / Velocity [m/s]
2 1.5 1 0.5 0 −0.5 −1 −1.5 −2 −2.5 0
10
20
30
40
50
Time [s]
Fig. 3.5: Phase relation between excitation force and actual velocity with the optimal control law
the generation power trajectories are almost the same. It turns out that the range here is only between -23 kW and 58 kW, so significantly less energy storage is required. However, the average generated power is less than with the optimal control law, namely 17.45 kW for “MPC Power” and 17.42 for “MPC Traj”. But that is still a large value given that the constraints restrict the allowed ranges for position, velocity and generator force by more than half. All in all it can be stated that both MPC formulations respect the constraints as desired and that the power formulation yields almost the same results as the trajectory formulation for the monochromatic wave data. The wave data modelling the incoming waves is then modified. The new input signal is a superposition of three sinewaves with amplitudes between 0.1 m and 0.3 m and period times between 1/2π s and 1/3π. Fig. 3.7 shows the simulation results under the same conditions apart from the wave data change. Again, it can be seen that all constraints are satisfied and that both formulations yield almost the same results. There is one interesting time period between 33 s and 43 s. In that period, the results of both formulations differ significantly. As can be seen on the bottom left, the velocity trajectory of “MPC Power” does not follow the commanded one, even though it is within the constraints. In this case, the power formulation calculates a completely different solution. Nevertheless, the constraints are satisfied and the average power is about the same. The average generated power for “MPC Traj” is 8.03 kW and for “MPC Power” is 7.87 kW . Interestingly, the optimal generated power with the optimal power law is only 3.72 kW. The wave signal in this simulation is a superposition of sine waves with different frequencies and amplitudes, and certainly not a monochromatic wave. Also, the signal
28
CHAPTER 3. ONE-BODY WEC Position of Wave and Buoy
Fgen (T) Fgen (P)
0.5
0.5
Force [N]
Position [m]
1
z (T) z (P) η
1
Generator Force
5
x 10
0 −0.5
0
−0.5
−1 −1 0
10
20
30
40
50
0
10
Time [s]
20
30
40
50
Time [s]
Buoy and commanded Velocity
Power Generation 80
z˙ (T) z˙ (P) z˙ com
1 0 −1 −2 0
Pgen (T) Pgen (P)
60
Power [kW]
Velocity [m/s]
2
40 20 0 −20
10
20
30
40
50
0
10
Time [s]
20
30
40
50
Time [s]
Fig. 3.6: Comparison between “MPC Traj” (T) and “MPC Power” (P) in case of monochromatic sinusodial wave as input signal.
is chosen so that there is no clear dominant wave frequency, making it one of the worst cases for this control law. This is the reason why the optimal power law is not optimal in this case. It is very interesting to see that the MPC algorithm of “MPC Traj” finds almost the same optimum trajectory for Fgen as “MPC Power”, even though “MPC Traj” tries to follow the “wrong” optimal reference trajectory. In summary, two different observations can be made: First, the formulation by means of the optimal trajectory yields good results even if the optimal trajetory is by far not the optimal one. That is an intersting result which was not expected in the first place. Secondly, the formulation by means of the direct power optimization leads to almost the same results as the other formulation. The advantage here is that no reference trajectory is necessary. Therefore, this approach can be easily applied to the linear two-body case, since an optimal velocity can not be calculated to easily fit the two-body case yet.
29
3.5. RESULTS
Position of Wave and Buoy 1.5
0.5 0 −0.5
0
−0.5
−1 0
Fgen (T) Fgen (P)
0.5
Force [N]
Position [m]
1
z (T) z (P) η
1
Generator Force
5
x 10
−1 10
20
30
40
50
0
10
Time [s]
Buoy and commanded Velocity
40
50
Power Generation 60
z˙ (T) z˙ (P) z˙ com
Pgen (T) Pgen (P)
40
Power [kW]
2
Velocity [m/s]
30
Time [s]
3
1 0 −1 −2 0
20
20 0 −20
10
20
30
40
50
−40 0
Time [s]
10
20
30
40
50
Time [s]
Fig. 3.7: Comparison between “MPC Traj” (T) and “MPC Power” (P) in case of a superposition of sine waves as input signal.
30
CHAPTER 3. ONE-BODY WEC
Chapter 4 Linear Model Predictive Control of Two-Body WEC This chapter presents a hydrodynamic two-body model of the point absorber in the time domain, which emphasizes the special mooring configuration and the formulation as a state space model. The next section deals with the model predictive control algorithm and extents the work done in the previous chapter. Additionally, two different approaches to handle possible infeasibility problems are shown before giving an insight into the implementation of the two-body MPC. Finally, simulation results of both approaches with irregular wave data are compared. The influence of horizon and step time of the control algorithm is also discussed by means of simulation results. Finally the influence of the mooring parameter for power generation is shown before the robustness behaviour is determined.
4.1
Two-Body Model
In this chapter, the point absorber is considered as a two-body device. This assumption is more realistic and better reflects the actual physical device. The left side of Fig. 4.1 illustrates the L10 point absorber developed at the Oregon State University. The right side shows a schematic two-body diagram which is used to derive the model. The floating body is called float or buoy, and the damping body is a combination of a spar and a ballast tank (in the following just called spar). Furthermore, the spar is moored to the sea floor in order to dampen the spar’s motion. Also, the two-body model is restricted to heave motion.
4.1.1
Equations of Motion
Much work has already been done on modelling two-body point absorbers. The herein presented modelling follows the work of Eidsmoen [13] and Ruehl [37]. In order to derive
32
CHAPTER 4. TWO-BODY WEC z Float PTO
Spar
x Tank
mooring
Sea Floor
Fig. 4.1: L10 Wave Energy Coverter [38] (left) and schematic diagram (right).
the equations of motion for the two bodies in the time-domain, some assumptions are made: • The formulation is based on Linear Wave Theorie (LTW). • Frequency dependent parameters of the L10 are assumed to be constant. • The spar is taut-moored with three cables. The pretension of the cables is induced
by the buoyancy forces of the spar. Therefore, the spar’s buoyancy force can be neglected.
• The radiation forces are assumed to be linear and no convolution terms are used to calculate it.
As can be seen on the right side of Fig. 4.1, x and z denote the position of the spar, respectively, of the buoy. Based on Newton’s law, the dynamic equations of the buoy can be written as Fgen + Fe1 − Fr1 − Fh1 − Fr12 = m1 z¨,
(4.1)
where z¨(t) is the buoy acceleration, m1 is the mass of the buoy, Fe1 is the excitaton force induced by the incoming waves, and Fgen is , as in the one-body model, the force produced by the power take-off system. Similar to the one-body case, the forces are Fr1 = A1 z¨(t) + b1 z(t), ˙
(4.2)
2 Fh1 = gρπrbuoy z(t) = k1 z(t),
(4.3)
where A1 is the buoy’s added mass, b1 is its viscous damping and k1 is the hydrostatic stiffness. Furthermore, there is a coupled radiation force resulting from the interaction of the two bodies. Here, the force Fr12 is acting from the spar on the buoy and can be written as Fr12 = A12 x¨.
(4.4)
33
4.1. TWO-BODY MODEL The same procedure can be applied to the second body. It follows that −Fgen + Fe2 − Fr2 − Fr21 − Fm = m2 x¨ Fr2 = A2 x¨ + b2 x˙ Fr21 = A21 z¨ Fm = 3K(sin α)2 x = Km x.
(4.5) (4.6) (4.7) (4.8)
There is no restoring or buoyancy force, since it is assumed that this force causes the pretension in the mooring cables. The mooring force is dependent on the adopted configuration. In this thesis, a taut-moored buoy with 3 cables is considered as shown in Fig. 4.2. The three mooring
Fig. 4.2: Mooring configuration with three cables (left) and schematic diagram for force derivation (right).
cables are configured like a tripod and there is an angle α with respect to the sea floor. According to Hook’s law, the mooring force in the direction of the cable can be written as F = −K(L′ − L),
(4.9)
where K is the cable stiffness. The vertical mooring force is then Fm = −3F sin(α).
(4.10)
Using the cosine theorem yields r π Fm = −3K L2 + x2 − 2Lx cos( + α) − L sin(α). (4.11) 4 However, this is a nonlinear law for the mooring. It is expected that linearization around x = 0 yields good results, since x 0. The same result with the opposite sign follows for x < 0. Even Fm,lin = Fm (x = 0) +
for large displacements such as 0.5 m, the difference between the nonlinear and the linear mooring force calculation is only 0.04% for a cable length of 170 m. Hence, the linear case is a very good approximation and the mooring force can be written as in (4.8).
4.1.2
State Space Form and Discretization
In order to control the two-body point absorber via MPC, a discrete time state space representation is necessary as well. Here, the formulation as a state space model is not straightforward. The problem lies in the coupling terms within the coupled radiation calculation. This coupling occurs in the second derivative, thus a state space model can not be formed without reformulating first. Equations (4.1)-(4.8) can be rewritten as Fgen + Fe1 − b1 z˙ − k1 z − A12 x¨ = z¨ (4.13) m1 + A1 −Fgen + Fe2 − b2 x˙ − Km x − A21 z¨ = x¨. (4.14) m2 + A2 Transforming these equations in state space form requires plugging (4.13) into (4.14) and vice versa to get rid of the coupling terms A21 z¨ and A12 x¨. The equations can be restated as A12 A21 = Fgen + Fe1 − b1 z˙ − k1 z z¨ m1 + A1 − m2 + A2 {z } | me1
+
A12 A12 b2 A12 Km A12 Fgen − F e2 + x˙ + x m 2 + A2 m2 + A2 m2 + A2 m2 + A2 (4.15)
A21 A21 A12 = −Fgen + Fe2 − b2 x˙ − x¨ m2 + A2 − Fgen m1 + A1 m1 + A1 {z } | me2
(4.16)
A21 b1 A21 Km A21 F e1 + z˙ + z − Km x m1 + A1 m1 + A1 m1 + A1 Now, it is possible to find a state space formulation with the state vector h iT h iT and the initial conditions x0 = 0 0 0 0 . Moreover, it is x = z z˙ x x˙ −
35
4.2. MPC FORMULATION
assumed that the entire state is measurable. The system in state space form can be written as x˙ = Ax + BFgen + B1 Fe1 + B2 Fe2 ,
(4.17)
y = x, where
A=
1
− mke1
− mbe1
1
A21 k1 me2 (m1 +A1 )
1 me1
0 − me
A21
2 (m1 +A1 )
A12 Km me1 (m2 +A2 )
0
A21 b1 me2 (m1 +A1 )
Km −m e
+
0
2
− mbe2
2
A12 me1 (m2 +A2 )
− m1e −
A12 b2 me1 (m2 +A2 ) 1
0 1 me1
0
0
0
0 1
0
B=
B1 =
0
A21 me2 (m1 +A1 )
2
,
(4.19)
0
A12 − me (m 2 +A2 ) 1 B2 = . 0
,
(4.18)
(4.20)
1 me2
In what follows, the inputs Fgen , Fe1 and Fe2 are denoted by u, v and w, respectively. With the same discretization approach as in Chapter 2, the following discrete-time system can be obtained x0 ∈ R4 ,
xk+1 = Ad xk + Bd uk + Bd1 vk + Bd2 wk , y k = xk ,
(4.21)
where Ah
Ad = e , Z h eAφ dφ B1 , Bd 1 =
Bd = Bd2 =
0
4.2
Z
Z
h
eAφ dφ B,
(4.22)
eAφ dφ B2 .
(4.23)
0 h 0
MPC Formulation
The MPC formulation for the two-body case is very similar to the direct power formulation in the one-body case. There is no information about an optimal reference trajectory, so this formulation is also the appropriate way to deal with the two-body problem. The optimization problem for the two-body case can be defined as min J(xk , uk ),
xk ,uk
(4.24)
36
CHAPTER 4. TWO-BODY WEC
where J(xk , uk ) = −
N X k=1
subject to
−q z˙k − x˙ k uk−1 −r u2k−1 , {z } |
(4.25)
−Pgen
xk+1 = Ad xk + Bd uk + Bd1 vk + Bd2 wk
(4.26)
|zk − xk | ≤ cp , ∀k = 1, . . . , N
(4.27)
|z˙k − x˙ k | ≤ cv , ∀k = 1, . . . , N
(4.28)
|uk | ≤ cu , ∀k = 0, . . . , N − 1,
(4.29)
One main difference is that the position and velocity constraints are now referred to the relative position and relative velocity between buoy and spar, since the stroke length is determined by the relative displacement. Also, the generated power is calculated by the product of the generator force and the relative velocity. Again, the following augmented vectors are used h iT h iT X = xT1 , xT2 , ..., xTN , U = u0 , u1 , ..., uN −1 , h iT h iT V = v0 , v1 , ..., vN −1 , W = w0 , w1 , ..., wN −1 .
(4.30) (4.31)
Solving system (4.26) and substituting in itself yields
X = Jx x0 + Ju U + Jv V + Jw W, where
Ad
2 Ad Jx = Ad 3 , . .. N Ad
Jv(w)
and
Bd
0
(4.32)
0
Ad B d Bd 0 A 2B Ad B d Bd Ju = d d . .. Ad N −1 Bd Ad N −2 Bd
Bd1(2) 0 0 Ad Bd1(2) Bd1(2) 0 A 2B Ad Bd1(2) Bd1(2) = d d1(2) .. . N −1 Ad Bd1(2) Ad N −2 Bd1(2) dim(Ju , Jv , Jw ) = (4N × N ), dim(Jx ) = (4N × 4).
...
... ... ...
... ... ... ... ...
0
0 0 , .. . Bd 0 0 0 , .. .
(4.33)
(4.34)
. . . Bd1(2)
(4.35)
37
4.3. INFEASIBILITY HANDLING Now, the objective function (4.25) can be reformulated in vector form as T ˆU J(X , U) = q X(2) − X(4) U + U T R T ˆU = q S1 X − S2 X U + U T R T
= qX (S1 |
T
(4.36)
ˆU − S )U + U R {z 2 } T
T
S
ˆ = diag(r) with dim(R) ˆ = (N × N ) and the matrices S1 , S2 with the matrix R
with dim(S1 , S2 ) = (N × 4N ). S1 and S2 extract the second state (buoy velocity),
respectively, fourth state (spar velocity) from the augmented state vector X . 0 1 0 0 0 0 ... 0 0 0 0 0 0 0 0 1 . . . 0 0 0 S1 = . .. ... , . . .
(4.37)
0 0 0 0 0 0 ... 1 0 0
0 0 0 1 0 0 0
... 0
0 0 0 0 0 0 0 1 . . . S2 = ... .. .
0 .. . .
(4.38)
0 0 0 0 0 0 0 0 ... 1
Using (4.32), the objective function can be formulated as a quadratic function that depends on U, namely,
ˆ U + q xT0 Jx T + V T Jv T + W T Jw T S U. J(U) = U T qJu T S + R
(4.39)
The constraint formulation as well needs just a few adjustments based on the one-body case. Provided that all constraints can be written as in (3.44) and (3.45), the MPC can be written as min J(U),
(4.40)
U
where J(U) = U T qJu T S + R U + q xT0 Jx T + V T Jv T + W T Jw T S U,
(4.41)
subject to
"
DD
#
"
d
#
U≤ . (4.42) E D Ju e − E D Jx x 0 − E D Jv V − E D Jw W This problem is also a convex optimization problem, as shown in Appendix A.
4.3
Infeasibility handling
In general, constrained optimization problems makes searching for the optimal solution more difficult. Furthermore, a serious problem can arise: the optimization problem is infeasible. That means there exists no solution within the boundaries. If this situation occurs, the optimization stops without a solution and the control loop is open. This
38
CHAPTER 4. TWO-BODY WEC
breakdown of the control can cause serious damage to the device and should be avoided. There are different possibilities to handle infeasibility problems. Most approaches try to deal with the infeasibility problem after it has occured. A simple method consists of setting the control signal equal to a constant value. Also, another type of control law, for example a P, PI or PID controller can be used in the event that the MPC does not find a solution. However, both implementations can not guarantee stability of the system within the constraints, since the system was already in a critical state when the breakdown occured. Another method prioritizes constraints and does not consider the least important ones while solving the problem again. However, it can not be assured how much the constraints get violated. This thesis introduces two approaches regarding hard and soft constraints. The first approach tries to find a feasible solution by relaxing some constraints after an infeasibility condition occurs, and resolves the optimization problem for the expanded region. Here the generator force limit cu and the velocity constraint cv can be relaxed, since overspeed and over-rated generator conditions can be tolerated for the short term [7]. The constraints for the new optimization can be defined by |zk − xk | ≤ cp ,
∀k = 1, . . . , N
(4.43)
|z˙k − x˙ k | ≤ rv cv , ∀k = 1, . . . , Nr
(4.44)
|z˙k − x˙ k | ≤ cv ,
∀k = Nr + 1, . . . , N
(4.45)
|uk | ≤ ru cu , ∀k = 0, . . . , Nr − 1
(4.46)
|uk | ≤ cu ,
(4.47)
∀k = Nr , . . . , N
(4.48) where rv and ru are the relaxation factors for velocity and generator force, respectively, and Nr denotes the number of time steps into the future until the non-relaxed constraints are valid again. This helps improve the stability attributes of the MPC, since it allows to violate the limits only at the beginning of the horizon. The second approach tries to prevent the occurance of infeasibility problems from the start. The original hard constraints are replaced by soft constraints through the introduction of positive slack variables ǫp , ǫv . In order to penalize violations of constraints, quadratic terms including the slack variables are inserted into the objective function that yields the new optimization problem min Jnew (U, ǫp , ǫv ),
(4.49)
Jnew (U, ǫp , ǫv ) = J(U) + wp ǫ2p + wv ǫ2v
(4.50)
U ,ǫp ,ǫv
where
39
4.4. IMPLEMENTATION subject to |zk − xk | − ǫp ≤ cp ,
∀k = 1, . . . , N
(4.51)
|z˙k − x˙ k | − ǫv ≤ cv , ∀k = 1, . . . , N
(4.52)
|uk | ≤ cu , ∀k = 0, . . . , N − 1 0 ≤ ǫ p ≤ c ǫp ,
(4.53)
0 ≤ ǫ v ≤ c ǫv .
(4.54)
The constraints cǫp , cǫv are hard upper bounds for the slack variables which are neccessary, in particular, for the position. The weighting factors wp , wv can be arbitrarily choosen and weight the penalty for violating the position and the velocity constraints, respectively.
4.4
Implementation
The simulation is also implemented using Matlab/SIMULINK. A diagram of the simulation model is shown in Fig. 4.3. Here, different types of irregular wave data obtained from [36] can be loaded. The prediction for the incoming wave elevation is calculated in the “Prediction calculation” block, where it is assumed that the prediction is perfect. In the next block, the excitation forces for both bodies are calculated by means of causal and noncausal impulse response functions from [36] as well. In that work, the impulse responses are obtained using frequency domain hydrodynamic analysis with ANYSY/AQUA [1] for the L10 wave geometry and parameters. The actual MPC is implemented similarly to the one-body case in a Level-2 Matlab S-Function, and obtains the future predictions for the excitation forces of both bodies and the current state of the device as inputs. The actual optimization problem is solved by means of the Matlab QP-solver “quadprog”. In order to check the status of the optimization, several other outputs are given to find, for instance, the points where the infeasibility catching is active. Moreover, it is assumed that all states are measurable.
eta
Fe1
Prediction _2
12:34 Simulation Time
eta1
Prediction calculation
UY
F_e_1
UY
F_e_2
F_gen_opt
Fe2
Excitation Force Determination
pos_buoy
Exitflag Optimization
F_gen TwoBody_diff F_gen_opt Fgen
z
pos_buoy
zdot
vel_buoy
vel_buoy
−20000 PTO
vel_buoy
Iteratations first optimization pos_plate
F_e_1 vel_plate
Fe_1 x
pos_plate
xdot
vel_plate
vel_diff Add
vel_plate F_e_2
Fe_2
Exitflag Infeasibility catching MPC
RealPlant
Fig. 4.3: SIMULINK simulation model for two-body WEC control with MPC.
40
CHAPTER 4. TWO-BODY WEC For the system parameters of the L10, the results in [36] obtained by ANSYS AQWA
[1] analysis are used. An overview is shown in Table 4.1, where the mooring constant Variable
Value
m1
2625.3 kg
m2
2650.4 kg
A11
8866.7 kg
A22
361.99 kg
A12
361.99 kg
A21
361.99 kg
b1
5,000 N/(m/s)
b2
50,000 N/(m/s)
k1
96,743 N/m
Km
300,000 N/m
Tab. 4.1: System parameter of the two-body L10 WEC.
which is within the is based on an angle α = 60◦ and a cable stiffness K = 133,333 N m common range for the stiffness of mooring cables, according to [10].
4.5
Results
Computer simulations are used to validate the MPC algorithm and its performance regarding the generated power. Furthermore, the formulation with relaxed constraints is compared to the soft constraints formulation. Another focus is on the behaviour when changing the parameters ∆t and Thor . In the end, the MPC is tested with several changes of parameters only in the “RealPlant” block, in order to simulate system uncertainties. A frequency spectrum from January 2009 of the NDBC Offshore buoy 46050 which is deployed off the coast of Oregon west of Newport [26] is used to generate a time series for the wave elevation. In all simulations this time series is used as the input signal to calculate the excitation forces. Six different cases are treated in this section. First, the controller behaviour of “Case 1” and “Case 2” are compared. The parameters for “Case 1” are listed in Table 4.2. This MPC is implemented with the relaxed constraint formulation. This means that the generator force constraint of 80,000 N can be violated in case an infeasibiliy condition occurs up to a value of 120, 000 N and the velocity constraint up to 2
m . s
In “Case
2”, the MPC is implemented by means of soft constraints. This formulation has no parameter values for Nr , rv and ru . The changed and the new parameters are listed in Table 4.3. To get comparable parameters, the hard constraint for the position is also set to 1 m consisting of 0.9 m soft constraint and the slack variable constraint of 0.1 m,
41
4.5. RESULTS Variable Explanation
Value
∆t
sample time for optimization
0.1 s
Thor
Optimization horizon
3s
N
Number of values
30
T
Simulation Time
100 s
q
Power weighting factor
1
r
Input weighting factor
0
cp
Relative position constraint
1m
cv
Relative velocity constraint
1 m/s
cu
Generator force constraint
80,000 N
Nr
Relaxation number
10
rv
Relaxation factor velocity
2
ru
Relaxation factor generator force
1.5
Tab. 4.2: MPC parameter values for the formulation with relaxed constraints (Case 1). Variable
Explanation
Value
cp
Relative position constraint
0.9 m
c ǫp c ǫv
Slack variable position constraint Slack variable velocity constraint
0.1 m 1 m/s
wp wv
Weighting factor position slack Weighting factor velocity slack
108 106
Tab. 4.3: MPC parameter values for the formulation with soft constraints (Case 2).
and the same holds for the velocity constraint. Furthermore, the weighting factor for the position slack is much greater than for the velocity, in order to penalize violations of the position constraint harder than violations of the velocity constraint. To get a first insight into the motion of the incoming waves and of the device, Fig. 4.4 shows wave, buoy and spar position for “Case 1”. The wave elevation η has a significant wave height of 4 m. It can be seen that the buoy position slightly lags the wave elevation and can not achieve the same amplitude as the wave elevation due to the constraints and the induced generator force. In general the spar position is within ±0.3 m.
More interesting is how often an infeasibility problem occurs, and Fig. 4.5 shows
that the optimization problem is infeasible at several points. A consecutive time span within which the optimization problem is not feasible is around 48 s. This is an interesting point to look at in the next figure as well. Fig. 4.6 shows the comparison between “Case 1” and “Case 2”. On the top left it can be seen that the hard position constraint is never violated in either formulation. The trajectory in “Case 2” tries to respect the soft constraint, since the penalty for violations is strong. The opposite can
42
CHAPTER 4. TWO-BODY WEC 2.5 z x η
2 1.5
Position [m]
1 0.5 0 −0.5 −1 −1.5 −2 −2.5 0
20
40
60
80
100
Time [s]
Fig. 4.4: Wave elevation, buoy position and spar position for the relaxed constraints formulation (Case 1). The wave data is from January 2009 of the NDBC Offshore buoy 46050.
be seen for the relative velocity on the bottom left. Since the penalty for violations is less, “Case 2” tends to violate this soft constraint more often. The generator force constraint in “Case 2” is never violated which is the purpose of the formulation. Only in the relaxed formulation is this constraint allowed to relax and it seems to be used at several instants in time, in particular around 48 s. At that time, the genarator force almost reaches its hard limit. On the bottom right, the power generation is shown and there is no noticeable difference between both lines in the figure. However, the average power of “Case 1” is 26.13 kW and is only 24.62 kW for “Case 2”. A possible reason for this is, for example, the possibility of having a larger generator force in “Case 1”. Moreover, in “Case 2” the MPC tries to respect the soft constraint of 0.9 m and thus limits the usable stroke displacement which leads to less power extraction. That can be reduced by lowering the weighting factor for the position slack. It should also be noted that the optimization problem in “Case 2” is always feasible, which is a big advantage in regards to the computational effort of the MPC. Indeed, the soft constraints formulation contains two more optimization variables (ǫp , ǫv ) which slightly increases the computation time, but in many cases it is possible to operate the MPC without infeasibility problems. Since a MPC is supposed to run in real time, a recalculation of the optimization problem as in “Case 1” is not always possible. So it is adventageous to have a formulation which is in general able to calculate the solution without recalculating as in the soft constraints formulation. However, it has to be said that even the soft constraints formulation can yield infeasibility problems (though
43
4.5. RESULTS
Infeasibility Catch
1
−2
0
20
40
60
80
100
Time [s]
Fig. 4.5: Frequency of the infeasibility problem for Case 1. 1 means original problem is feasible. −2 means original problem is infeasible.
rarely), since the slack variables are also subject to constraints. So, for the worst case scenario it makes sense to combine both formulations so that it is possible to relax the generator force constraint, even in “Case 2” if an infeasibility problem occurs. Next, two cases are compared in terms of the number of optimization variables. In both cases the soft constraints formulation of “Case 2” is used. In “Case 3” the horizon time is set to 6 s from 3 s and in “Case 4” the step time is doubled to 0.2 s. Therefore, the weighting factors for the slack variables are set to half of the original value in “Case 3” and are doubled in “Case 4” in order to get comparable conditions. It is remarkable that the results of “Case 3” are almost the same as in “Case 2”. Likewise, the average generated power of 24.44 kW is similar as well. Therefore, the results of both cases are compared directly, without comparison with “Case 2” which can be seen in Fig. 4.7. Here as well, no noticeable difference can be seen and the average power in “Case 4” is 24.34 kW. The advantage of the larger time step is that the number of optimization variables decreases. “Case 2” contains 32, “Case 3” 62 and “Case 4” only 17 optimization variables. This makes a huge difference regarding the computation time, since the computational effort increases with more optimization variables. “Case 4” yields almost same results, even though the sample time is doubled. But one important difference can be seen in Fig. 4.7, especially in the relative velocity plot. The trajectory calculted with the larger time step tends to oscillate much more which can cause problems (particulary in case of large inaccuracies between model and the actual device). However, it can be seen that a horizon of 3 s is sufficient for this type of wave data.
44
CHAPTER 4. TWO-BODY WEC
Position Difference of Buoy and Plate Case 1 Case 2
0.5
0.5
0
0
−0.5
−0.5
−1
−1
0
Case 1 Case 2
1
Force [N]
Position [m]
1
Generator Force
5
x 10
20
40
60
80
100
0
20
Time [s]
Velocity Difference of Buoy and Plate
80
100
Power Generation 120
Case 1 Case 2
1.5
Case 1 Case 2
100 80
1
Power [kW]
Velocity [m/s]
60
Time [s]
2
0.5 0 −0.5
60 40 20 0
−1 −1.5 0
40
−20 20
40
60
80
−40 0
100
20
Time [s]
40
60
80
100
Time [s]
Fig. 4.6: Comparison between relaxed (Case 1) and soft (Case 2) constraints formulation of the MPC.
Position Difference of Buoy and Plate
Velocity Difference of Buoy and Plate 2
Case 3 Case 4
1
1.5 1
Velocity [m/s]
Position [m]
0.5
0
−0.5
0.5 0 −0.5 −1
−1 0
Case 3 Case 4
20
40
60
80
Time [s]
100
−1.5 0
20
40
60
80
100
Time [s]
Fig. 4.7: Comparison between the soft constraints formulation with the changed parameter Thor = 6s (Case 3) and the same formulation with changed parameter ∆t = 0.2s (Case 4).
45
4.5. RESULTS
Next, only the influence of the mooring parameter Km on the system peformance is discussed. Therefore, the mooring constant is changed from Km = 300,000 N/m (“Case 2”) to Km = 50,000 N/m (“Case 5”). Both cases are simulated with the soft constraints formulation. The results are shown in Fig. 4.8. The position and velocity trajectories Position Difference of Buoy and Plate Case 2 Case 5
0.5 0 −0.5 −1 0
Case 2 Case 5
1
Force [N]
Position [m]
1
Generator Force
5
x 10
0.5 0 −0.5
20
40
60
80
−1 0
100
20
Time [s]
Velocity Difference of Buoy and Plate
80
100
Power Generation 120
Case 2 Case 5
1.5
Case 2 Case 5
100 80
1
Power [kW]
Velocity [m/s]
60
Time [s]
2
0.5 0 −0.5
60 40 20 0
−1 −1.5 0
40
−20 20
40
60
Time [s]
80
100
−40 0
20
40
60
80
100
Time [s]
Fig. 4.8: Comparison between the soft constraints formulation with mooring constant Km = 300,000 N/m (Case 2) and the same formulation with Km = 50,000 N/m (Case 5).
on the left do not differ significantly. However, a large difference can be seen in the generator force. The generator force in “Case 5” is within its limits as well. But it is in general far away from the limits and much smaller than in “Case 2”. Due to the smaller generator force, it can be seen, that the generated power is also much smaller, namely, only 11.91 kW in comparison to 24.62 kW on average. It stands out as well that the required reactive power is less which is advantageous regarding the required energy storage. Fig. 4.9 shows the spar position for both simulations. Because of the extreme change of the mooring constant, the mooring force decreases significantly and thus, the spar motion is within a larger range up to ±1.3 m in comparison to only ±0.3
m in “Case 2”. This is the reason why the relative position and velocity between buoy and spar are naturally less with smaller mooring forces. So, the generator force does not need to exceed its limits to respect position and velocity constraints. Hence, the generated power decreases, since the relative motion difference decreases. Generally, it can be stated that a stronger mooring increases the power generation. However, it
46
CHAPTER 4. TWO-BODY WEC 1.5 Case 2 Case 5
Position Spar [m]
1
0.5
0
−0.5
−1
−1.5 0
20
40
60
80
100
Time [s]
Fig. 4.9: Comparison of Spar position between the soft constraints formulation with mooring constant Km = 300,000 N/m (Case 2) and the same formulation with Km = 50,000 N/m (Case 5).
has to be pointed out that a stronger mooring brings the generator faster to its limits. Therefore, especially for rough wave conditions, it might be advantageous to have a looser mooring to prevent violations of constraints and physical damage of the device. Moreover, in case there is a chance to alter the mooring configuration during operation, that can be used to improve the power generation in mild wave climates by increasing the mooring constant and in rough conditions, the mooring constant could be reduced to protect the device from damage. Finally, the MPC is tested regarding inaccuracies between model and actual device. Hence the parameters of the system model inside the MPC and the excitation force calculation remain the same, whereas the parameters of the actual system change. This definitely does not test the behaviour of every uncertainty. It is only a check regarding parameter uncertainties and does not include unmodelled dynamic uncertainties or non-linearities in the actual device. Nevertheless, it is a good first try to get an idea of the robustness attributes of the proposed MPC. The changed real plant parameters for “Case 6” and the unchanged model parameters are displayed in Table 4.4. Basically, the mooring and damping forces of the assumed actual WEC are reduced and the added mass and restoring force of the buoy are augmented. In order to compare the results with the unchanged model, “Case 2” and “Case 6” are shown together in Fig. 4.10, where there are significant differences in the velocity and position trajectories. But it is remarkable that the constraints are not (or only slightly) violated. However the average generated power is only 18.59 kW. Overall it can be stated that the MPC with the changed system parameters is still
47
4.5. RESULTS Variable
Real Plant
Model
A11
9,500 kg
8866.7 kg
A22
320 kg
361.99 kg
A12
380 kg
361.99 kg
A21
340 kg
361.99 kg
b1
4,000 N/(m/s)
5,000 N/(m/s)
b2
42,000 N/(m/s)
50,000 N/(m/s)
k1
120,000 N/m
96,743 N/m
Km
210,000 N/m
300,000 N/m
Tab. 4.4: Comparison of real plant and model parameters. Position Difference of Buoy and Plate Case 2 Case 6
0.5
0.5
0
0
−0.5
−0.5
−1
−1
0
Case 2 Case 6
1
Force [N]
Position [m]
1
Generator Force
5
x 10
20
40
60
80
100
0
20
Time [s]
Velocity Difference of Buoy and Plate
80
100
Power Generation 120
Case 2 Case 6
1.5
Case 2 Case 6
100 80
1
Power [kW]
Velocity [m/s]
60
Time [s]
2
0.5 0 −0.5
60 40 20 0
−1 −1.5 0
40
−20 20
40
60
Time [s]
80
100
−40 0
20
40
60
80
100
Time [s]
Fig. 4.10: Comparison between the original soft constraint formulation with original WEC parameters of Table 4.1 (Case 2) and the same formulation with the changed WEC parameters of Table 4.4 for the L10 (Case 6).
stable and remarkably robust with respect to the system constraints. For a better overview, all discussed cases with explanation are shown in Table 4.5.
48
CHAPTER 4. TWO-BODY WEC
Case
Explanation
1
Relaxed constraint formulation
2
Soft constraint formulation
3
Soft constraint formulation with larger step time (∆T = 0.2 s)
4
Soft constraint formulation with larger horizon (Thor = 6 s)
5
Soft constraint formulation with smaller mooring constant (Km = 50,000 N/m)
6
Changed system parameter for robustness check Tab. 4.5: Overview Cases.
Chapter 5 Nonlinear Model Predictive Control of a Nonlinear Two-Body WEC Model The goal of this chapter is to propose a nonlinear model predictive controller to deal with a nonlinear WEC model. It presents first a brief introduction to nonlinear model predictive control (NMPC). Following, the nonlinear model in [36] with convolution terms and a nonlinear mooring law is discussed and compared to the model without convolution terms. Also, the discretization of this nonlinear system and the implementation of the NMPC is a part of this chapter. By means of simulation results of controlling the nonlinear model, the proposed NMPC is validated and compared to the linear MPC from Chapter 4.
5.1
Introduction to Nonlinear Model Predictive Control
Linear MPCs are well known and have been applied since the 1970’s, whereas the attention to NMPC arises only in the 1990’s [16]. Linear MPC theory is quite mature today and system theoretic attributes as stability and optimality are well addressed [22]. Also, many different industrial MPC applications can be found. The case is different with NMPC. While theoretical characteristics are well discussed, industrial applications are difficult to find [5]. Linear MPCs and NMPCs have basically the same concepts. In general, the NMPC problem is formulated as solving a finite horizon optimal control problem which is subject to the constraints and to the system dynamics [16]. The basic principle is the same as described in Chapter 2 and can also be seen in Fig. 2.1. Whereas only linear models are used to predict future system behaviour in linear MPC techniques,
50
CHAPTER 5. NONLINEAR MPC
in NMPC, nonlinear models are used. Furthermore, even system constraints can be considered nonlinear. The motivation for NMPC arises, since there are systems which are inherently nonlinear and can not be described adequately by a linear model. Also, control and process requirements have increased. In these cases, the use of linear MPC is inadequate and NMPC is a promising alternative to deal with inherently nonlinear systems. A fundamental problem of NMPC schemes is that the constrained optimization problem needs to be solved within a specified time limit. In the case of linear MPC, the problem is convex and for the class of LQP, proven optimization algorithms exist to solve the problem efficiently. NMPC requires the solution of a nonlinear problem, though. In general, these problems are non-convex, thus it can not be assured to find the global optimum. Additionally, the solution can be computationally expensive. Therefore, it is important to exploit the special structure of each problem to obtain a real-time feasible optimization problem. This work does not focus on real-time applicability. The focus is on the qualitative performance of NMPC regarding a nonlinear WEC model whose nonlinearity results from a nonlinear mooring force. So, this work attempts to establish if NMPC is advantageous to use for controlling the nonlinear WEC compared to the linear MPC.
5.2
Nonlinear Two-Body Model
The proposed nonlinear two-body model follows the work in [36]. In addition to the proposed two-body model in Chapter 4, impulse response functions are used in order to calculate the radiation forces. These functions are obtained from [36]. For further information regarding the impulse response functions, the readers are referred to [36] as well. Using impulse response functions yields convolution terms in the expressions for the radiation forces. Furthermore, a highly nonlinear mooring force law, [36, 21], is used. The equation of motions for the two bodies are: Fgen + Fe1 − Fr1 − Fh1 − Fr12 = m1 z¨
(5.1)
−Fgen + Fe2 − Fr2 − Fr21 − Fm = m2 x¨
(5.2)
with the following terms for the forces Fr1 = A1 z¨(t) +
Z
t −∞
fr1 (t − τ )z(τ ˙ ) dτ,
2 Fh1 = gρπrbuoy z(t) = k1 z(t), Z t fr12 (t − τ )x(τ ˙ ) dτ, Fr12 = A12 x¨(t) + −∞
(5.3) (5.4) (5.5)
51
5.2. NONLINEAR TWO-BODY MODEL Z t Fr2 = A2 x¨(t) + fr2 (t − τ )x(τ ˙ ) dτ, −∞ Z t fr21 (t − τ )z(τ ˙ ) dτ, Fr21 = A21 z¨(t) +
(5.6) (5.7)
−∞
Lnl . (5.8) L2nl + x(t)2 The impulse response functions of the different radiation forces are denoted by fr . The Fm = 8Knl x(t) 1 − p
mooring system is based on the experimental mooring configuration in [21]. There, the buoy is moored with 8 cables to a static reference around the buoy. No mooring to the sea floor is assumed. Knl is the stiffness of one cable, and Lnl is the horizontal length from the reference to the buoy. The system expressed by (5.1)-(5.8) only differs from the proposed linear two-body system in Chapter 4 in the use of the convolution terms and in the mooring force. So, it is interesting to see how large the influence of the convolution terms is. Neglecting these terms and assuming the mooring force (5.8) yields the nonlinear model x˙ = Anl x + l(x) + Bu + B1 v + B2 w
(5.9)
where
Anl
=
0
1
0
− mke1
− mbe1
0
1
1
0
A21 k1 me2 (m1 +A1 )
0
0
A21 b1 me2 (m1 +A1 )
0
0
A12 b2 me1 (m2 +A2 ) 1
− mbe2 2
0 8A12 Knl me (m2 +A2 ) 1 − √ L2nl 2 x Lnl +x 1 l(x) = , 0 Lnl nl √ − 8K x 1 − 2 me 2 2
,
(5.10)
(5.11)
Lnl +x
and B, B1 and B are the same as in the linear two-body model from Chapter 4. Now there are two different models to compare. The first model (5.1)-(5.8) with convolution terms is called extended model, and the second model (5.9) is called reduced model. Fig. 5.1 shows the relative position trajectories for simulations with both models. This simulation uses wave data from buoy 40650 and applies a control law proportional to the velocity difference. The control law was chosen simply, because only the influence of the convolution terms should be determined. No difference between both results can be seen on the left side, whereas zooming in shows a slight difference on the right side which is only about 0.4%, even in the worst cases. Hence, it is reasonable to neglect the convolution terms and using the reduced model for the NMPC approach. In case there would be a larger influence of these terms, the convolution terms can be approximated by a linear state space model by methods described in [19, 39]. However, this is not necessary here.
52
CHAPTER 5. NONLINEAR MPC 2
1.772 extended reduced
extended reduced 1.771
1
Relative position [m]
Relative position [m]
1.5
0.5 0 −0.5 −1
1.769
1.768
1.767
−1.5 −2 0
1.77
20
40
60
80
100
1.766 45.8
45.82
45.84
Time [s]
45.86
45.88
45.9
Time [s]
Fig. 5.1: Comparison between the extended and the reduced model. The wave data is from January 2009 of the NDBC Offshore buoy 46050 [26]. Proportional control is applied. On both sides, the relative position is shown. On the right side only a small detail can be seen.
5.3
Discretization
In the following, a finite parameterization of the controls and constraints is used to find a direct solution of the optimization problem. Therefore, the system is supposed to be described as a discrete-time nonlinear state space model in the form xk+1 = f k (xk , uk , sk )
(5.12)
where the function f k (.) maps the current state xk , the manipulable input uk and the unmanipulable inputs sk to the next state xk+1 . According to [20], a nonlinear system x˙ = f (x(t)) + g(x(t))u(t)
(5.13)
with real analytic vector fields f (x(t)) and g(x(t)) can be discretized by an approximate sampled-data representation under Zero-Order Hold assumption by xk+1 = xk +
M X l=1
where A[l+1] (x, u) =
A[l] (xk , uk )
hl , l!
δA[l] (x, u) (f (x) + g(x)u) δx
(5.14)
(5.15)
with l = 1, 2, 3, . . .. Here M denotes the order of the discretization, and h is the sample time. In simulation results it turns out that a discretization of order 1 (comparable to Euler forward method) is not appropriate for the proposed NMPC approach. In fact, it yields unreasonable results. Due to this fact, the discretization order is chosen to be M = 2 in the
5.4. PROBLEM FORMULATION AND IMPLEMENTATION
53
ensuing work. With the nonlinear system (5.9), it follows that A[1] = Anl x + l(x) + Bu + B1 v + B2 w,
(5.16)
A[2] = [Anl + l′ (x)] · [Anl x + l(x) + Bu + B1 v + B2 w].
(5.17)
Using (5.14), after some manipulations the discrete-time system can be described by h2 ′ h2 l (xk )Anl xk + l′ (xk )l(xk ) 2 2 ′ +Bd u + Bd1 v + Bd2 w + l (xk )(Bu + B1 v + B2 w),
xk+1 = Anl,d1 xk + Anl,d2 l(xk ) +
(5.18)
where h2 h2 2 Anl,d1 = I + hAnl + Anl , Anl,d2 = hI + Anl , 2 2 h2 h2 Bd = hB + Anl B, Bd1 = hB1 + Anl B1 , 2 2 h2 Bd2 = hB2 + Anl B2 , 2 with the jacobian matrix 0 0 0 0 L3nl 12 Knl 0 0 m 8A(m 0 · 1 − (L2 +x2 )1.5 e1 2 +A2 ) ′ nl k,(3) l (xk ) = . 0 0 0 0 3 L nl nl 0 0 − 8K 0 · 1 − 2 2 1.5 me ) (L +x 2
5.4
nl
(5.19) (5.20) (5.21)
(5.22)
k,(3)
Problem Formulation and Implementation
The NMPC is also implemented using Matlab/SIMULINK. The simulation model has the same structure as the linear MPC for the linear two-body model. The problem solver “quadprog” is not appropriate for this case, since this solver can only handle LQPs. Therefore, the solver “fmincon” of Matlab is used which can handle nonlinear problems with nonlinear constraints as well. In the following, the solution of the nonlinear optimization problem using Matlab is outlined. First, the optimization problem can be stated as min
xk ,uk ,ǫ1 ,ǫ2
J(xk , uk , ǫ1 , ǫ2 ),
(5.23)
where J(xk , uk , ǫ1 , ǫ2 ) =
N X k=1
subject to
q xk,(2) − xk,(4) uk−1 + r u2k−1 + wp ǫ21 + wv ǫ22 , xk+1 = f k (xk , uk , vk , wk )
(5.24)
(5.25)
54
CHAPTER 5. NONLINEAR MPC |xk,(1) − xk,(3) | − ǫ1 ≤ cp ,
∀k = 1, . . . , N
(5.26)
|xk,(2) − xk,(4) | − ǫ2 ≤ cv , ∀k = 1, . . . , N
(5.27)
|uk | ≤ cu , ∀k = 0, . . . , N − 1 0 ≤ ǫ p ≤ c ǫp ,
0 ≤ ǫ v ≤ c ǫv ,
(5.28) (5.29)
with the same denotation as in Chapter 4. The problem is implemented using the following optimization vector iT h p = xT1 u0 xT2 u1 . . . xTN uN −1 ǫ1 ǫ2 ,
(5.30)
with 5N + 2 variables and it is straightforward to express the objective function by means of p. The solver “fmincon” can handle box constraints, linear and nonlinear inequality and equality constraints. The slack variables (5.29) and the input constraints (5.28) are considered as box constraints which yield 2N + 4 equations. The position (5.26) and velocity (5.27) limits are considered as linear inequality constraints which yield 4N equations. Additionally, the system dynamics (5.25) (or rather (5.18)) need to be included as constraints, here, as nonlinear equation constraints with 4N equations. In the following, Thor = 3s and ∆T = 0.1s, and thus the horizon length is 30 steps. In summary, the optimization problem consists of 152 optimization variables and 304 equations for the constraints. This is a large problem for online solving of the problem within a step time of ∆T = 0.1s. However, as stated above, the goal of this work is to determine the qualitive behaviour of the system controlled by the NMPC in comparison with a linear MPC and not to discuss real-time applicability.
5.5
Results
In the next simulations the WEC parameters from Table 4.1 are used, and a cable stiffness of Knl = 840,000 N/m and a cable length to the reference of Lnl = 1.7 m for the mooring force calculation is assumed. The system which is controlled is always the nonlinear system (5.1)-(5.8) with convolution terms. Two different controllers are proposed to control this system: The NMPC as described in Section 5.4 and a linear MPC which is implemented as shown in Chapter 4. The NMPC includes the nonlinear model without convolution terms as prediction model, whereas the linear MPC is based on the linear model with linear moring force constant Km . Both controllers are supposed to control the nonlinear system with nonlinear mooring force and with convolution terms. The results with the NMPC and the linear MPC are denoted by NL and LIN, respectively. The simulations and control parameters are shown in Table 5.1. In what follows, the NMPC is validated and compared to the linear MPC with different values
55
5.5. RESULTS Variable
Explanation
Value
∆t
sample time for optimization
0.1 s
Thor
Optimization horizon
3s
N
Number of values
30
T
Simulation Time
50 s
q
Power weighting factor
1
r
Input weighting factor
0
cp
Relative position constraint
1m
cv
Relative velocity constraint
0.9 m/s
cu
Generator force constraint
50,000 N
c ǫp
Slack variable position constraint 0.1 m
c ǫv
Slack variable velocity constraint
1 m/s
wp
Weighting factor position slack
108
wv
Weighting factor velocity slack
106
Tab. 5.1: NMPC/MPC parameter values.
of Km through simulation. Fig. 5.2 shows the nonlinear and linear mooring forces for certain values of Km for spar displacements between ±0.6 m. It can be seen that the nonlinear law is highly
nonlinear and a linearization would not work. Still, there are many possibilities to choose Km . In the following, three different cases are regarded. One (Km = 50,000 N/m) which underestimates the nonlinear mooring, another (Km = 350,000 N/m) which assumes a large mooring force and a third choice of Km so that the linear MPC generates the maximal power for the current wave data. In the beginning, a time-series from the NDBC Umpqua buoy 46229 which is deployed off the coast of Oregon north of Reedsport [26] is used as the wave data. In Fig. 5.3 the comparison between the NMPC and the linear MPC with Km = 100,000 N/m is shown. First, it can be noted that the proposed NMPC yields reasonable results, the constraints are satisfied and a large amount of power is generated. For this choice of Km , the linear MPC yields its best result. Obvious differences are between 15 s and 30 s. However, these differences in the trajectories are not significant and have no effect on the generated power, since the NMPC yields a average power of 9.45 kW and the linear MPC 9.41 kW. More significant differences can be seen in Fig. 5.4, where the results with the NMPC are compared to two different cases. The linear MPC which assumes a high mooring force yields very different trajectories and worse results than the two cases. Due to the assumption of high mooring forces, especially for small displacements, the MPC attempts to prevent the bouy’s large displacements, in order to satisfy the constraints.
56
CHAPTER 5. NONLINEAR MPC 5
2.5
x 10
NL K=50000 K=125000 K=350000
2
Mooring Force [N]
1.5 1 0.5 0 −0.5 −1 −1.5 −2 −2.5 −0.6
−0.4
−0.2
0
0.2
0.4
0.6
Position of Spar [m]
Fig. 5.2: Nonlinear and linear mooring forces for different choices of Km .
Since the actual mooring force is smaller, the relative motion between spar and buoy is reduced. Thus, the average generated power for this case is only 5.42 kW. The MPC with a smaller mooring constant yields clearly better results. The average generated power is 8.39 kW which is less than in the optimal case. However, this choice has also certain advantages. As can be seen, the position and generator force constraints are very well satisfied, where the trajectories remain well within the constraint bounds. Also, the required reactive power is much less. Regardless of the values of Km the simulation with NMPC produces the best power generation for this wave data. In the following, the controller performance with different wave data, where the generator force limit is also increased to 80,000 N, is tested. The wave data is from the NDBC Offshore buoy 46050 with a little rougher wave climate. Even though the dominant frequency is smaller, the wave amplitudes are significantly larger. The optimal value of Km for this wave data is 150,000 N/m. Fig. 5.5 shows the comparison between NMPC and the linear MPC for the new wave data, where some differences can be seen. The NMPC trajectories run far away from their constraints, thus, the entire power potential is not used. Hence, the average generated power for the NMPC is only 17.98 kW, compared to 22.8 kW for the linear MPC. Even though the parameter Km is perfectly chosen and there exist many values which do not yield such results, it is remarkable that the results with the linear MPC are significantly better. It can thus be concluded that the NMPC has probably certain weak points. One reason might be the discretization method of only order of 2, whereas the linear discretization is exact. Furthermore, in some cases the nonlinear solver might only find the local optimum which yields non-optimal results. In future research, these problems need to be given attention. Nevertheless, the results with NMPC satisfy the constraints and
57
5.5. RESULTS Position Difference of Buoy and Plate 1.5
NL Km=100000
6
Force [N]
Position [m]
8
NL Km=100000
1
Generator Force
4
x 10
0.5 0
4 2 0 −2
−0.5
−4 −1 0
10
20
30
40
−6 0
50
10
20
30
Time [s]
Time [s]
Velocity Difference of Buoy and Plate
Power Generation
40
50
2 NL K =100000
80
m
0.5 0
20 0
−1
−20 10
20
30
40
m
40
−0.5
−1.5 0
NL K =100000
60
1
Power [kW]
Velocity [m/s]
1.5
50
−40 0
Time [s]
10
20
30
40
50
Time [s]
Fig. 5.3: Comparison between the simulation results with the NMPC and the linear MPC with Km = 100,000 N/m. The wave data is from the NDBC Umpqua buoy 46229.
are reasonable, even though they are not optimal. It is interesting to see, for which values of Km the linear MPC yields the best results and how sensitive this controller is with respect to changes of Km . In Fig. 5.6, the generated power depending on the value of Km for both wave datas is shown. As stated above, the result with the NMPC is the best for the milder wave climate. However, with a mooring constant of about Km = 100,000 N/m the results with the linear MPC are almost the same. Especially for higher mooring constants, the average generated power decreases significantly. For the rougher wave climate, there is a large range of values of Km , where the results with the linear MPC are even better than with the LMPC. There is a considerable range of values for Km in both cases, where changes of Km do not affect the generated power very much. These regions are near the optimum value. It can be observed that the optimum performance for the rougher wave climate is achieved by a larger value of Km , namely, 150,000 N/m. This can be explained by the spar position. The spar position for the milder wave climate is around ±0.35 m,
whereas it is ±0.42 m for the rougher climate, since the wave exerts more force to
the bodies. Due to the larger displacements, the mooring force increases significantly because of the highly nonlinear mooring force as seen in Fig. 5.2. In this plot it can
58
CHAPTER 5. NONLINEAR MPC Position Difference of Buoy and Plate 1.5
NL Km=50000
6
K =350000 m
K =350000 m
4
Force [N]
Position [m]
8
NL Km=50000
1
Generator Force
4
x 10
0.5 0
2 0 −2
−0.5 −1 0
−4 −6 10
20
30
40
50
0
10
20
30
Time [s]
Time [s]
Velocity Difference of Buoy and Plate
Power Generation
40
50
2 80
NL K =50000
1.5
m
0.5 0 −0.5
m
K =350000 m
40 20 0 −20
−1
−40
−1.5 0
NL K =50000
60
K =350000
1
Power [kW]
Velocity [m/s]
m
10
20
30
40
50
Time [s]
0
10
20
30
40
50
Time [s]
Fig. 5.4: Comparison between the simulation results with the NMPC and the linear MPC with Km = 50,000 N/m and Km = 350,000 N/m. The wave data is from the NDBC Umpqua buoy 46229.
be seen that with the optimal value of Km the linear mooring force is similar to a best-fit line for the nonlinear mooring force. In this case, higher nonlinear mooring forces occur, thus the best-fit line is naturally a line with a larger slope, which results in a larger mooring constant Km .
59
5.5. RESULTS
Position Difference of Buoy and Plate 1.5
NL K =150000
m
Force [N]
Position [m]
NL K =150000
10
m
1 0.5 0 −0.5 −1 0
Generator Force
4
x 10
5
0
−5
10
20
30
40
50
0
10
20
30
Time [s]
Time [s]
Velocity Difference of Buoy and Plate
Power Generation
40
50
2
1
Power [kW]
Velocity [m/s]
100
NL Km=150000
1.5
0.5 0 −0.5
NL Km=150000
50
0
−1 −1.5 0
10
20
30
40
−50 0
50
Time [s]
10
20
30
40
50
Time [s]
Fig. 5.5: Comparison between the simulation results with the NMPC and the linear MPC with Km = 150,000 N/m. The wave data is from the NDBC Umpqua buoy 46050.
15
25 LIN(K )
LIN(K )
12.5
m
NL
Average Generated Power [kW]
Average Generated Power [kW]
m
10
7.5
5
2.5
0 0
22.5
NL
20 17.5 15 12.5 10 7.5
0.5 1 1.5 2 2.5 3 3.5 5 Assumed Mooring Stiffness K [N/10 ]
4
5 0
0.5 1 1.5 2 2.5 3 3.5 5 Assumed Mooring Stiffness K [N/10 ]
4
Fig. 5.6: Generated Power depending on Km compared to the generated power with the NMPC. On the left side with wave data from buoy 46050 and on the right side with data from buoy 46229.
60
CHAPTER 5. NONLINEAR MPC
Chapter 6 Conclusions and Future Work This thesis presented several different model predictive control approaches for controlling wave energy converters. It exclusively discussed point absorbers, though the general concepts could be applied to other types of WECs as well. In all types of control approaches the goal was to maximize the generated power while satisfying generator force as well as position and velocity constraints for the buoy. Also, different state space models and a nonlinear model for the point absorber has been derived. Two different MPC formulations were proposed to control a WEC modelled as a one-body device. Both formulations have been compared with each other and to a well-known optimal control law. It was successfully demonstrated that it is possible to control the one-body WEC efficiently without using an optimal velocity trajectory for the buoy velocity. Moreover, both formulations yielded almost the same results. Since the calculation of an optimal velocity trajectory is not always possible, especially in the two-body case, this was an important achievement in order to contol a two-body device. Thus, it could be shown that a two-body model can be controlled by a MPC scheme as well. This fact has been demonstrated by simulation results with real wave data from an offshore buoy. Additionally, the infeasibility problem of optimization problems was discussed and two different approaches to handle this problem were proposed, namely, by relaxing the constraints and by introducing soft constraints. Both approaches yielded good results, where the generated power with the proposed soft constraint formulation was in general a little less than with the relaxed constraint formulation. However, it could be shown that the soft constraint formulation handles the infeasibility problem better and avoids its occurance from the start which is also computationally cheaper. Also, it has been shown that the proposed horizon and step time is a good choice to deal with the applied wave data and that the MPC is quite robust against parameter changes of the WEC. Additionally, the influence of the mooring of the spar on the system performance
62
CHAPTER 6. CONCLUSIONS AND FUTURE WORK
was under examination. The results showed that a better mooring of the spar leads in general to a higher power generation. However, the danger that the constraints can not be maintained increases with better mooring, since the relative motion between buoy and spar increases. Thus, the mooring configuration has an essential influence on the WEC performance. So the idea arises that variable mooring, which means adjusting the mooring stiffness during operation, could be used to improve the system’s safety and performance, if variable mooring is technically possible. Many assumptions for the linear models has been made that are not always appropriate. That is why a nonlinear MPC was also proposed in this work to deal with mooring and radiation force nonlinearities. It has been demonstated that the NMPC is able to keep the point absorber states within its limits. In order to evaluate the performance, the results were compared to a linear MPC with varying mooring constant. It was shown that the performance with the linear MPC is even better under certain wave conditions and for certain values of the mooring constant. It has been concluded that the performance of the NMPC could be improved by a better discretization and a specialized optimization method. Overall, it can be stated that MPC is a very promising and adequate approach to control point absorber wave energy converters, since MPCs are able to maximize the generated power while taking physical constraints into account, which is exactly the challenge for controlling WECs. However, much research still remains be done in the future to bring ocean wave energy applications to commercialication and expanded implementation. Therefore, this work suggests several areas for future research. In all simulations, a perfect prediction of the incoming wave elevation is assumed and hence a perfect prediction of the excitation forces. The implementation of prediction algorithms such as a Kalman filter or auto-regression methods can enhance this work and can identify the influence of prediction errors. Additionally, this will impact optimal energy capture performance. Furthermore, it was assumed that the entire state of the device is measured which is not always possible. If not, an observer has to be implemented to identify all states, since the MPC needs the full state information. The proposed two-body model is still not validated against experimental data and all parameters arise only from AQUA analysis. In case the model structure is not appropriate to describe the point absorber well enough, more precice hydrodynamics including frequency dependence or other non-linearities have to be introduced. The proposed NMPC sets a framework for dealing with these nonlinearities. However, some improvements regarding discretization and solver need to be considered in future research. Additionally, the field of real-time applicability needs to get center stage for all
63 types of proposed MPCs. The real-time implementation of the linear MPCs with two, respectively, four states and 30 instants of time does not cause large difficulties regarding online computation time. More work regarding solver and suitable implementation platform needs to be done in the nonlinear case.
64
CHAPTER 6. CONCLUSIONS AND FUTURE WORK
Appendix A Explantion: Direct Power Formulation is a Convex Problem Optimization problems can have many different structures. The problem (4.41,4.42) which needs to be solved in the Two-Body MPC algorithm is a QLP. For this class of optimization problems there exist well known and suitable methods to solve it [27], like the interior-point or active-set method. However, finding the global optimum requires convex problems. Thus, the local optimum is the global optimum. The usual tracking formulation of MPC is inherently convex. In the following, it is explained why the power formulation as well yields a convex problem. For a general QLP with objective function 1 (A.1) J = xT Qx + f T x, 2 if Q is a positive definite matrix then the optimization problem is convex. Moreover, if all eigenvalues of Q are positive, then Q is positive definite. For problem (4.41,4.42) Q = qJu T S + R. Since q > 0 and R is a diagonal matrix with elements r ≥ 0, it
is sufficient to consider only Ju T S. This matrix is an upper triangular matrix and therefore, it is positive definite if all values on the diagonal are positive. Bd (2) − Bd (4) ⋆ ... ⋆ .. (2) (4) B − B . ⋆ d d Ju T S = . . .. ..
∅
(A.2)
Bd (2) − Bd (4)
Since the optimization problem is convex if Bd (2) − Bd (4) > 0, where (2), (4) denote the second and forth component of the vector Bd , respectively. The continuous state
ii
APPENDIX A. CONVEX PROBLEM
space vector B has always the form 0 + B= 0. −
After discretization,
Bd =
∞ X i=1
Therefore, Bd has at least the form
Ai−1 B
(A.3)
hi . i!
(A.4)
⋆ + Bd (A.5) ⋆, − even if the discretization is not exact, since the first component Bh of the sum is the dominant factor. It follows that Bd (2) − Bd (4) is always positive and hence the opti-
mization problem is convex and can be solved, for example, using Matlab’s quadprog solver which deals with convex quadratic problems.
Bibliography [1] ANSYS Inc. ANSYS AQWA 13.0. 275 Technology Drive, Canonsburg, PA. [2] C. K. Alexander and M. N. O. Sadiku. Fundamentals of Electric Circuits. Tata McGraw-Hill Publishing Company Limited, 2007. [3] E. Arnold. Numerische Methoden der Optimierung und Optimalen Steuerung. Institute for Systemdynamics, University of Stuttgart, 2010. script. [4] A. Babarit and A. Cl´ement. Optimal latching control of a wave energy device in regular and irregular waves. Applied Ocean Research, 28(2):77 – 91, 2006. [5] T. Badgwell and S. Qin. Nonlinear predictive control - theory and practice. The Institution of Electrical Engineers, 2001. [6] T. Brekken, A. von Jouanne, and H. Y. Han. Ocean wave energy overview and research at oregon state university. In Power Electronics and Machines in Wind Applications, 2009. PEMWA 2009. IEEE, pages 1 –7, jun. 2009. [7] T. K. Brekken. On model predictive control for a point absorber wave energy converter. IEEE. [8] J. Buijs, J. Ludlage, W. V. Brempt, and B. D. Moor. Quadratic programming in model predictive control for large scale systems. Technical report, University Leuven, Department of Electrical Engineering, 2002. [9] E. Camacho and C. Bordons. Model Predictive Control. Springer, 2003. [10] S. K. Chakrabarti. Handbook of Offshore Engineering. Elsevier, 2005. [11] C. Chen. Linear System Theorie and Design. Oxford University Press, 1999. [12] A. Cl´ement et al. Wave energy in europe: current status and perspectives. Renewable and Sustainable Energy Reviews, 6(5):405 – 431, 2002. [13] H. Eidsmoen. Simulation of a slack-moored heaving-buoy wave-energy converter with phase control. PhD thesis, Norwegian University of Science and Technology, Trondheim, Norway, 1996.
iv
BIBLIOGRAPHY
[14] A. F. O. Falc˜ao. Wave energy utilization: A review of the technologies. Renewable and Sustainable Energy Reviews, 14(3):899 – 918, 2010. [15] J. Falnes. Ocean Waves and Oscillating Systems, Linear Interaction Including Wave-Energy Extraction. Cambridge University Press, 2002. [16] R. Findeisen and F. Allg¨ower. An introduction to nonlinear predictive control. In 21st Benelux Meeting on Systems and Control, Veldhoven, 2002. [17] Global
Greenhouse
Warming.
Wave
Dragon.
http://www.
global-greenhouse-warming.com/wave-dragon.html, sep. 2011. [18] J. Hals, T. Bjarte-Larsson, and J. Falnes. Optimum reactive control and control by latching of a wave-absorbing semisubmerged heaving sphere. ASME Conference Proceedings, 2002(36142):415–423, 2002. [19] E. Jefferys. Simulation of wave power devices. Applied Ocean Research, 6(1):31 – 39, 1984. [20] N. Kazantzis, K. T. Chong, J. H. Park, and A. G. Parlos. Control-relevant discretization of nonlinear systems with time-delay using taylor-lie series. Journal of Dynamic Systems, Measurement, and Control, 127(1):153–159, 2005. [21] Y. Li, Y. Yu, M. Previsic, E. Nelson, and R. Thresher. Numerical and experimental investigation of a foating point absorber wave energy converter under extreme wave condition. [22] D. Mayne, J. Rawlings, C. Rao, and P. Scokaert. Constrained model predictive control: Stability and optimality. Automatica, 36(6):789–814, 2000. [23] S. McArthur and T. Brekken. Ocean wave power data generation for grid integration studies. In Power and Energy Society General Meeting, 2010 IEEE, pages 1 –6, jul. 2010. [24] C. Moler and C. V. Loan. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. Siam Review, 45:3–49, 2003. [25] A. Muetze and J. Vining. Ocean wave energy conversion - a survey. In Industry Applications Conference, 2006. 41st IAS Annual Meeting. Conference Record of the 2006 IEEE, volume 3, pages 1410 –1417, oct. 2006. [26] National Data Buoy Center. http://www.ndbc.noaa.gov/, sep. 2011. [27] J. Nocedal and S. J. Wright. Numerical Optimization. Springer Verlag, 2006.
v
BIBLIOGRAPHY
[28] Northwest National Marine Renewable Energy Center (NNMREC). The Limpet OWC by Wavegen. http://nnmrec.oregonstate.edu/wavegen-owc, sep. 2011. [29] Ocean Power Technologies (OPT). Global Resources, sep. 2011. [30] Pelamis Wave Power. http://www.pelamiswave.com, sep. 2011. [31] H. Polinder and M. Scuotto. Wave energy converters and their impact on power systems. In Future Power Systems, 2005 International Conference on, pages 9 pp. –9, nov. 2005. [32] T. Pontes et al. The european wave energy resource. In Proceedings of the 3rd European Wave Energy Conference, Patras, Greece, 1998. [33] J. Prudell, M. Stoddard, E. Amon, T. Brekken, and A. von Jouanne. A permanentmagnet tubular linear generator for ocean wave energy conversion. Industry Applications, IEEE Transactions on, 46(6):2392 –2400, nov.-dec. 2010. [34] C. V. Rao, S. J. Wright, and J. B. Rawlings. Application of interior-point methods to model predictive control. Journal of Optimization Theory and Applications, 99:723–757, 1998. [35] M. Richter. Modell-praediktive Trajektorienplanung zur aktiven Seegangskompensation. Student thesis, Institute for Systemdynamics, 2010. [36] K. Ruehl. Time-domain modeling of heaving point absorber wave energy converters, including power take-off and mooring. Master’s thesis, Oregon State University, Mechanical Engineering, 2011. [37] K. Ruehl, T. Brekken, B. Bosma, and R. Paasch. Large-scale ocean wave energy plant modeling. In Innovative Technologies for an Efficient and Reliable Electricity Supply (CITRES), 2010 IEEE Conference on, pages 379 –386, sep. 2010. [38] E. Rusch. Catching a wave, powering an electrical grid? In Smithsonian Magazine, 2009. [39] R. Taghipour, T. Perez, and T. Moan. Hybrid frequency–time domain models for dynamic response analysis of marine structures. Ocean Engineering, 35(7):685 – 705, 2008. [40] U.S. Energy Information Administraion.
International Energy Outlook 2011.
http://www.eia.gov/, sep. 2011. [41] J. Vining and A. Muetze. Economic factors and incentives for ocean wave energy conversion. Industry Applications, IEEE Transactions on, 45(2):547 –554, mar.apr. 2009.
vi
BIBLIOGRAPHY
[42] Voith Hydro Wavegen Limited. Wavegen brochure. http://www.wavegen.co. uk/, sep. 2011. [43] Wave Dragon. http://www.wavedragon.net, sep. 2011. [44] A. Wills and W. Heath. EE03016 - interior-point methods for linear model predictive control. Technical report, University of Newcastle, Australia, 2003.
