Multi-Objective Predictive Control Optimization with Varying Term Objectives: A Wind Farm Case Study

Abstract

This paper introduces the incentive of an optimization strategy taking into account short-term and long-term cost objectives. The rationale underlying the methodology presented in this work is that the choice of the cost objectives and their time based interval affect the overall efficiency/cost balance of wide area control systems in general. The problem of cost effective optimization of system output is taken into account in a multi-objective predictive control formulation and applied on a windmill park case study. A strategy is proposed to enable selection of optimality criteria as a function of context conditions of system operating conditions. Long-term economic objectives are included and realistic simulations of a windmill park are performed. The results indicate the global optimal criterium is no longer feasible when long-term economic objectives are introduced. Instead, local sub-optimal solutions are likely to enable long-term energy efficiency in terms of balanced production of energy and costs for distribution and maintenance of a windmill park.

1. Introduction

When it comes to optimization strategies, advanced control methodologies such as model based predictive control (MPC) is of great industrial relevance [1,2,3]. For large scale systems, its variant as distributed MPC has great added value in terms of numerical optimization and computational efficiency [4,5]. Additionally, it requires a significantly lower amount of information than full multivariable MPC, hence the optimization can be accelerated. The cost function implemented in such MPC schemes is usually tailored upon the specific objectives of the process at hand. Less academic but highly relevant in practice objectives such as performance degradation, failure monitoring and nesting, implementation and training costs are making the MPC an even more appealing strategy for large scale processes. Supply chain optimization is a great example for successful stories of economic and environmental related objective based MPC applications [6,7].

Wide area control systems have been recently undergoing a revision of concepts and relevance going beyond the classical output performance [8,9]. Such examples of revisited control objectives can be found in large scale power systems [10], office buildings [11], traffic control [12], sustainability in company management [13], and distribution networks [9]. A generic feature is to search solutions for optimal operation from the decentralized to distributed control systems and combinations thereof. The features of tomorrow’s control systems are based on system-wide control, plug-and-play control, measurement driven, adaptive and reconfigurable architectures.

An application where environmental and economic objectives are core in the long-term cost management and return of investment policies is the area of renewable energy, e.g., wave energy [14]. Example of large scale interacting sub-systems with independent control but global optimization policies is the case of windmill parks, either land- or water-based. Due to their negative impact on aesthetics and noise near populated areas, land based windmill parks lose their popularity as the marine parks gather higher interest [15]. Satellite-based synthetic aperture radar measurements allows a most optimal location along with best geometrical placement of such parks for the function of wind speed and its energy content [16,17]. These issues are now mature and the next at hand objective for optimal operation is derived from evidence-based effects of wind power systems on marine life and related maintenance costs—the latter being obviously quite different from land-based parks [18]. Reports of various agencies and policy makers indicate that long-term and large-scale impacts have cumulative effects and a knowledge gap is identified.

Following the exponential growth of wind capacity since 1996 to now, sufficient data is now available to comprehensively evaluate their efficiency, by taking into account the effective energy production and related costs on the long-term. Recent studies evaluating the distribution of costs for wind power systems indicate that offshore based parks have about double costs with respect to land-based costs [19]. Of the total cost distribution, maintenance and safe operation at lower power production are estimated to be a staggering 30%–40%, and thus the major source of costs [18,20]. The same distribution of costs applies also to a levelised cost of electricity from wind power systems: operation and maintenance remain leaders in the cost indicators [21].

Maintenance costs are related to the operation of the wind parks, in terms of environmental conditions (saline air, saline water, etc) and safety related operational decisions (wind speed and stability of the construction). There are other effects such as variability of the wind speed, affecting the grid where the power is introduced, further enlarging the cost-related issues to the grid itself [22]. Choice of instrumentation is also important to avoid operational problems in the grid [23,24,25]. Weak grids will operate at a lower power than the capacity for power generation from the wind park itself, and the choice of control strategy (voltage- or power-based criterium) will greatly affect the impact [26]. Limitations in storage capacity will impose limits in peak voltage as a safety margin.

In this paper, we perform a feasibility study of a concept for multi-objective optimization (MOO) scheme integrated with MPC and applied on a relevant study case: a windmill park. The combined scheme addresses optimization criteria on a priority-based algorithm, resulting in a multi-objective optimization priority MPC (MOOP-MPC ). Economic costs are taken into account along with safety related limitations [27]. Comparison to a global priority optimization formulation is performed to analyse the best strategy for long-term impact. An existing model for a windmill park and updated cost function is used in Matlab/Simulink. Here, the estimation of wind speed is the basis for control feedback law and a real-time estimation scheme is proposed. The hypothesis tested in this study is whether the choice of the MOOP objectives will have an effect on the global performance output and efficiency of the system.

2. Materials and Methods

Model predictive control is an advanced control strategy well established over the past decades and commonly employed in the so-called money making industry, where economic costs are highly relevant [28,29,30]. For processes consisting of multiple interacting sub-systems, with highly coupled dynamics, multivariable MPC algorithms are demanding in terms of model data availability and computational costs. Instead, a much lighter version, ignoring interactions (if weak enough) or used with decoupling matrices (if interactions are strong), is that of decentralized MPC. As a trade-off solution, with interaction information communicated among the sub-systems, is the distributed MPC [4,5], also a well established MPC strategy.

In this paper, we use a distributed MPC algorithm presented in [4,5], with a tailored optimization scheme described hereafter. For nonlinear systems, in our version of MPC, linearization of the process model is not necessary, in the condition that the step response of system dynamics matrix is updated at every sampling time and the step input to obtain it has the amplitude in the region of the expected steady state values of the controller output (due to nonlinear dynamics, if large input is used, the information matrix no longer has information upon the specific operation point currently used). This is a special version of the MPC, especially suitable for data driven formulation; see comprehensive details in [31,32].

2.1. Multi Objective Optimization with Priorities

The new approach for wide-area control of systems is to have decentralized and distributed control layers throughout the global system output to be controlled [10,33]. Assuming the inner/lower layer of control works adequately, the upper level can deal with objectives other than absolute output, or combinations thereof with economic and environmental objectives.

At this point, we introduce the prioritized multi-objective optimization (MO) algorithm. This is a simplified approach compared to those proposed in the literature [34,35,36,37,38]. As with any process, the safety constraint is set as a hard constraint, given limit value intervals for all input–output variables. If this condition is not satisfied, a pre-set of (suboptimal) safety values are given to the process operation units. This step is implemented as proposed in [39]. If this condition is fulfilled, then the next priority is to meet the product specifications, i.e., performance tolerance error intervals and/or maintenance costs are evaluated. Finally, if this is also fulfilled, then the optimization minimizes the control effort, i.e., enters energy saving operation of the plant. The performance and other long-term costs are soft constraints, i.e., they are tailored to fit the objective at hand and not to minimize a specific cost goal. This allows a much faster computational convergence while process operation remains active within safety bounds. The sequential (prioritized) flowchart is iterated at every sampling period and the computational time within iteration is recorded.

2.2. Area-Wise Wind Speed Estimation

Bat algorithms and learning machine algorithms are popular estimation tools for predicting wind speed on location of the wind turbine [26,40]. By adjusting the wind turbine speed, the system can operate at optimal rotational speed while achieving the maximal power. However, the operation of the turbine is influenced by wind speed and direction in a complex context given the influence of neighbouring turbines in a windmill park. A good wind estimator is desired to reduce the uncertainty that will directly affect the performance of the proposed controller in a neighbouring area. Figure 1 depicts the position of the wind estimator within the loop.

Drones or in general, unmanned aerial vehicles (UAV) can be used to fly over the park and estimate wind speeds at determined locations in order to update the information for the individual operation of the turbines depending on the time-varying wind conditions. The proposed wind estimation is derived below for the coordinate frames of the UAV.

Solutions to estimate the wind using on-board embedded wind sensors, are given in [41,42]. A possible drawback is that the on-board wind sensors use valuable payload that could be otherwise used. From aerodynamics, a nonlinear dependence of the UAV follows, by the wind speeds $u_w$, $v_w$, $w_w$, while the disturbances (external forces $F_{Xaero}$, $F_{Yaero}$, $F_{Zaero}$ and moments $L_{aero}$, $M_{aero}$, $N_{aero}$) enter linearly in the drone equations. Hence, the problems of estimation of wind velocities and disturbances can be formulated in two context conditions: piecewise constant or slowly varying. To estimate linear parameters varying models for wind, the literature mainly reports two groups of solutions: with/without airspeed sensor information. For instance, Qu et al. [43] uses airspeed sensor, roll angle, and sideslip angle to estimate the wind speed for small fixed-wing UAV. Rhudy et al. [44] estimated the wind field based on four formulations, combining the data coming from the Pitot -static tube, the global position system (GPS), the inertial measurement unit (IMU) and the angle of attack and sideslip vanes. Otherwise, Xing et al. [45] estimated the shear wind vector at low altitude using IMU and GNSS modules. Pappu et al. [46] used a Kalman filter based gust identification technique for estimating wind gusts. In our work, it is supposed that the estimation algorithm can use IMU (accelerometer, gyroscope) sensors augmented with a motion tracking system and rotor’s rotational velocity sensors.

Quadrotor linear velocities $(u,v,w)$ and accelerations are provided by the on-board accelerometer, which measures directly

$${\tilde{\Delta }}_a(X)={\Delta }_a(X)+{ϵ}_a$$

(1)

where ${ϵ}_a$ is the bounded measurement noise and

$${\Delta }_a(X)=\left[\begin{array}{c}{\Delta }_{au}(X)\\ {\Delta }_{av}(X)\\ {\Delta }_{aw}(X)\end{array}\right]=\left[\begin{array}{c}\dot{u}\\ \dot{v}\\ \dot{w}\end{array}\right]+\left[\begin{array}{c}p\\ q\\ r\end{array}\right]X\left[\begin{array}{c}u\\ v\\ w\end{array}\right]-\left[\begin{array}{c}-gsin\theta \\ gcos\theta sin\varphi \\ gcos\theta cosphi\end{array}\right].$$

(2)

From the gyroscope, which measures the rotational velocity in body frame with respect to the earth, the other state coordinates are measured

$${\tilde{\Delta }}_g(X)={\Delta }_g(X)+{ϵ}_g$$

(3)

where ${ϵ}_g$ is the bounded measurement noise and

$${\tilde{\Delta }}_g(X)={\left[pqr\right]}^T.$$

(4)

The IMU sensor is augmented with ground based cameras, used to estimate $(u,v,w,\theta ,\varphi )$ coupled with the gyroscope and accelerometer, enabling drone observability with respect to the inertial frame. The dynamics equation can be summarized as follows:

$$\left[\begin{array}{c}{\Delta }_a(X)\\ {\Delta }_g(X)\end{array}\right]=f_0(X,U,\omega )+\mathsf{\Omega }(\omega )d_w$$

(5)

where $f_0$ is known and $\mathsf{\Omega }$ needs to be estimated.

Define the predicted acceleration

$${\tilde{\Delta }}_a=f_{0a}+{\mathsf{\Omega }}_a{\widehat{d}}_w$$

(6)

allowing introduction of the error between the measured and predicted state acceleration as

$$e_a={\tilde{\Delta }}_a-{\widehat{\Delta }}_a={\mathsf{\Omega }}_a(d_w-{\widehat{d}}_w)+{ϵ}_a.$$

(7)

From Rios et al. [47], we use the finite-time estimation algorithm

$${\widehat{d}}_w={\gamma }_a{\mathsf{\Omega }}_a^T{[e_a]}^{{\alpha }_a},0<{\alpha }_a<1,{\gamma }_a\gg 0.$$

(8)

The Lyapunov function can be selected

$$V=\frac{1}{2{\gamma }_a}{|d_w-{\widehat{d}}_w|}^2$$

(9)

with ${||}^2$ denoting the norm and for the simplified dynamics of a quadrotor, the matrix ${\mathsf{\Omega }}_a$ is invertible and the solution is feasible and stable [47,48].

An adaptive observer equation model can be derived at this point

$$\begin{array}{c}{\widehat{\mathsf{\Omega }}}_g=f_{0g}+{\mathsf{\Omega }}_g{\widehat{d}}_w+l_g\mathrm{sign}({\tilde{\Delta }}_g-{\widehat{\Delta }}_g)\\ {\widehat{d}}_w={\gamma }_g{\mathsf{\Omega }}_g^T({\tilde{\Delta }}_g-{\widehat{\Delta }}_g)\end{array}$$

(10)

with $l_g$ and $\gamma \gg 0$ tuning parameters. Consider here the Lyapunov function

$$V=\frac{1}{2}\left({[{\Delta }_g-{\widehat{\Delta }}_g]}^2+\frac{1}{{\gamma }_g}{[d_w-{\widehat{d}}_w]}^2\right)$$

(11)

which has the derivative in time

$$\dot{V}=-l_g|{\Delta }_g-{\widehat{\Delta }}_g|$$

(12)

and is bounded for all $t>0$, with the state estimator converging asymptotically to the origin [47]. The wind estimation error converges as well to the origin, due to the persistent excitation in ${\mathsf{\Omega }}_g$.

Some convergence bottlenecks may appear with this algorithm as the computations occur. To increase robustness of convergence, the following algorithm is proposed (recall operations are element-wise):

$$\begin{array}{c}{\widehat{\Delta }}_g=f_{0g}+{\mathsf{\Omega }}_g{\widehat{d}}_w+l_g^{{}^{\prime }}\left({\tilde{\Delta }}_g-{\widehat{\Delta }}_g\right)+\mathsf{\Phi }{\widehat{d}}_w\\ \mathsf{\Phi }=-l_g^{{}^{\prime }}\mathsf{\Phi }+{\mathsf{\Omega }}_g\\ {\widehat{d}}_w={\gamma }_g^{{}^{\prime }}{\mathsf{\Phi }}^T{|{\tilde{\Delta }}_g-{\widehat{\Delta }}_g|}^{{\alpha }_g^{{}^{\prime }}}\end{array}$$

(13)

where $l_g^{{}^{\prime }}>0$, ${\alpha }_g^{{}^{\prime }}\in (0,1)$ and ${\gamma }_g^{{}^{\prime }}\gg 0$ are tuning parameters. The auxiliary matrix in this augmented model has the same dimension of ${\mathsf{\Omega }}_g$ and is limited for bounded values of ${\mathsf{\Omega }}_g$ and positive values of $l_g^{{}^{\prime }}$.

In the presence of measurement noise and time-varying wind speed conditions, Equation (10) will freeze its estimates, while (8) and (13) will estimate in an aggressive, respectively conservative manner. The fusion of the equations may seem appealing under the form

$${\widehat{d}}_w^{fusion}(t)=\frac{{\sum }_{i=1}^2{e}^{-k_i{v}_i^2(t)}{\widehat{d}}_w^i}{{\sum }_{i=1}^2{e}^{-k_i{v}_i^2(t)}}$$

(14)

with $k_i$ positive, a tuning parameter for convergence speed, and errors $v_1(t)=e_a(t)$, respectively $v_2(t)={\tilde{\Delta }}_g(t)-{\widehat{\Delta }}_g(t)$. This fusion algorithm has the worst estimation error given by the maximum of the two estimation errors for the algorithms (8) and (13).

Important fluctuations in the wind speed such as rapid turbulence will result in voltage fluctuations. The above described algorithm can estimate the basal envelope of the speed, while a faster component can be estimated by filtered white noise. The mean and standard deviation of the wind speed are linearly related with a constant k found experimentally from the park site:

$${\widehat{\sigma }}_w=k·\overline{v}.$$

(15)

Using the shaping filter from Suvire et al. [49]

$$H(j\omega )=\frac{K}{{(1+j\omega T)}^{5/6}}$$

(16)

we can set the time constant of the filter as

$$T=\frac{L}{\overline{v}}$$

(17)

with L the turbulence length (e.g., hundreds of meters). The constant of the filter is set by the condition of coloured noise with unit standard deviation for the wind speed values:

$$K=\sqrt{\frac{2\pi ·T}{B(0.5,0.3)·T_s}}$$

(18)

with $T_s$ the sampling period and B denoting the Bessel function.

The global estimation of the wind speed will dictate the decisions in the MOOP-MPC optimization algorithm.

2.3. Windmill Park Simulator

A recent comprehensive literature review of wind farm operation using distributed predictive control illustrates that wind speed is a core variable of information in the optimization of the farm output in terms of generated power [50]. The wake effect is mainly in stream with the geometry of the park and has a decaying effect towards the end area of its direction. Integration of wind mill parks has been discussed comprehensively in a tutorial [10] and models for power control assumed from Ugalde-Loo et al. [51], with Matlab codes available as in Sadamoto et al. [52]. Networked delayed control is relevant for this system, but it is neglected in this study [8].

The simplest manner to model the wind park is to consider it a uniform turbine model. This is useful for global optimization objectives with global wind speed conditions. The other extreme is to consider individual models for each turbine. Wind mill parks can be as large as hundreds of units, hence a significant computational load makes such models restricted to specific analysis such as coherence, correlation, and similar properties of wind effects on energy production. However, area related (i.e., group, or local) effects on the optimization scheme may be investigated if groups of wind turbines are considered bundled into a model, and the wind park consists of several of such lumped models. We use the model presented in Suvire et al. [49] and available in Matlab/Simulink, with values for realistic wind speeds from Degraer et al. [53], but in a simplified version.

The schematic in Figure 2 depicts the concept of area-wise simulation of the windmill park and possible interaction between the areas coming from the wind direction.

The model of each sub-system denoted in Figure 2 by its corresponding number is a simplified one from the one used in Zhang et al. [40] and Suvire et al. [49], which is easily implementable in Matlab/Simulink. The output is the turbine’s generated power. The time constant of the low pass filter depends on the average wind speed and can be assumed constant of time-varying. The effect of the wind fluctuations at rated power operation is filtered by a damped second order transfer function. Figure 3 schematically depicts the system to be controlled. This denotes a single grouped area of windmills in the large park. The model has been fitted as a first hand least square optimization algorithm to determine the parameters: $T_b=11$, $K_p=2.8$, $T=12.5$ and $d=97$, the latter term being related to the distance covered by one sub-system.

3. Results

To verify the result between a global and an area-wise MOOP with distributed MPC strategy, simulations are performed in Matlab/Simulink R2017a. The values of the turbine model are the same as in Zhang et al. [40], and the values for the distributed MPC algorithm are set to 30 samples for the prediction horizon, one sample as the control horizon, and no time delay. The sampling period is 100 milliseconds. The output of the system is the power. The tuning of the MPC follows the rules of thumb given in Ionescu et al. [54].

The MOOP has three parameters sequentially activated as in Figure 4. The safety limits for the input are 0–100 Volts, rotor flux for turbine between 0–1 and angle speed for rotor tracking between 10–40 rad/s. Current is limited to maximum 5 A. These are taken as hard constraints in the global optimization case, with a 100% weight. The maximum power is normalized to 1 pu, reaching maximum output at 10 m/s and safety shut-down resulting in zero power output at 20 m/s. For the simulation purposes, the default values from SimPowerSystems/Simulink Toolbox within Matlab are used, but with wind speed from Degraer et al. [53] with an average varying between 10–20 m/s. The length of the park is assumed to be 300 m. The total output power of the farm is obtained by adding the power from each sub-system area for which model from Figure 3 has been identified from the simulator data.

The performance term is evaluated as quadratic error between the maximal power extracted from the estimated wind speed and actual power from the system, during global optimization at a 70% weight. Maintenance costs are long-term costs and during global optimization they have a weighting factor in the second objective as a 30% of the total 100% weight—the rest is given to the performance as a short-term objective will have more influence on the cost variability. Alternatively, the short-term and long-term costs can be also implemented as a function of the prediction horizon: shorter intervals will give a faster convergence with large fluctuations, whereas longer prediction intervals will provide a more basal variation in the long-term objectives. A ratio of 2:5 is proposed for the short vs. long-term prediction horizon.

When exploring the environmental impact information such as in Degraer et al. [53], this is a cyclic signal, which can be modeled as a slowly moving average signal with increasing offset (i.e., cost). The minimization therefore is in terms of average values bringing them piecewise towards zero. The impact of this term during the global optimization case is 10%.

Figure 5 depicts the various timelines of the optimization parameters and the wind speed profile for the simulation is given in Figure 6.

The safety criterion is tested as being of the foremost importance in the operation of the park. Figure 7 illustrates the simulation of the park with a moment for wind speed above 20 m/s and temporary low windmill operation, resulting in lower output power—once the wind speed is recovered below the maximum allowed, the operation is resumed.

We compare now the two situations in variable speed wind conditions. First, all sub-system areas of the park are globally optimized, with a solidary cost function to take into account all output variables equally. Distributed MPC (dMPC) is still valid in this case, but the objective function is limited to neighbouring areas of each individual group of windmills. Second, the optimization is done individually per group (per sub-system), taking into account interactions from neighbouring areas, but the optimization cost function is a MOOP with weighting factors among the various types of costs. For the same time interval of 20 min, with same wind speed conditions, the power output for the two strategies is depicted in Figure 8.

In case of constant average speed conditions, the comparison between the two strategies is illustrated in Figure 9.

The obtained performance is summarized in

Table 1. Normalized units in percent for variable wind speed conditions. Distributed MPC (dMPC), Global and Multi-objective optimization priority (MOOP) methods.

Method	Performance	Total Power Output	Cost
Global dMPC	88	85	91
MOOP dMPC	63	71	68

and

Table 2. Normalized units in percent for constant wind speed conditions.

Method	Performance	Total Power Output	Cost
Global dMPC	90	95	91
MOOP dMPC	87	91	88

below for the two wind speed situations.

4. Discussion

The idea of using UAVs as a flying sensor is borrowed from precision agriculture applications [55] and emergency medicine [56], but its integration in windmill park utilisation still leaves many open challenges. For instance, the effect of wind gusts and wake on the UAV in flying mode may be destabilizing. Efforts to model such effects are numerous [57,58]. Also, the need for real time measurements of wind speed may not always be justified, but this again fits with the use of UAVs, and this is perhaps less appealing for on-board turbine instrumentation [59]. The advantage of using UAVs for sensing wind speed may be in their versatility to approach different areas of the park at different altitudes and at different moments in time, while still having relatively less costs when compared to on-board instrumentation. The problem of in flight stability has been addressed in manifold applications of UAV and it is considered a rather mature control problem and is beyond the objective of this paper. By contrast, the use of UAVs can be eliminated, at the cost of determining via measurements the wind speed, direction and other useful features for maximizing throughput of the windmill park. In this study we have particularly chosen the use of UAVs, motivated by the manifold applications and versatility of their use in limited environments.

The use of the simplified model for wind turbine and the simplified wind farm scheme have the advantage of allowing the simulations to convey information on the methodology of optimization and control, while keeping to a minimum the effects of difficult dynamics from the system itself. The disadvantage is that the analysis is limited to a concept, an incentive for further development, and does not claim to be a comprehensive feasibility study of the wide area control system at hand.

It is not yet clear whether the increase in output power with maximizing global power extraction from the park is actually more productive than the MOOP optimization scheme. In fact, they have to be analysed in the perspective of their maintenance costs, which as indicated, are higher for the global optimization case than in MOOP. The relatively lower power extraction of the MOOP case may be well justified in long-term economic investigation.

From the simulations we performed, it appears that under constant wind speed conditions, there is no significant difference among the two strategies. This is somewhat expected, as the variability among the sub-systems remains similar for uniform distributions of wind energy in time. As the area-based dynamics are simplified to common dynamics, there is no source of inter-system variability propagated throughout the park. By contrast, variable speed conditions will induce different cost function convergences, as the data used for optimization contains persistent excitation and different weights affect the cost function optimizer. In this case, we observe significant differences among the two strategies.

There are manifold opportunities for control-related items to be further investigated. Moreover, the following questions may be applicable to any wide-area control system.

For instance, how would comparing the global optimization of maximizing output power against a MOOP with a short-term cost function and a long-term cost function? Whereas by long-term cost function, one could use an economic cost of maintenance over several months, and impact on environment cost over several years.

Would it still be interesting/appealing to use the maximizing power strategy?

Or, perhaps, a better solution could be that such a strategy could be used alternatively with the MOOP short–long-term cost objectives?

These are merely a shortlisted enumeration of possible research directions. As mentioned in [10], we are witnessing a transition to a new infrastructure in a manifold of power system networks, creating a shift in energy supply from a centralized to a distributed network of energy supply and demand. Such changes are visible in other than the power sector as well. Climate changes and stringent environmental regulations affect the evaluations of cost-driven investments and this is visible in long-term rather than short-term criteria. The concept and incentive study presented in this paper suggest that operators are well motivated to be inclined to explore new control methods that go far beyond current system management.

5. Conclusions

This paper introduces an incentive for wide-area control systems, with a conceptual study for optimization of output power as a function of short-term and long-term impact objectives. A relevant windmill park case study is presented, in which the initial simulations indicate opportunities and challenges for both control and economic driven studies.