An extremal problem with applications to renewable energy production

Math. Model. Nat. Phenom. 16 (2021) 37 Mathematical Modelling of Natural Phenomena https://doi.org/10.1051/mmnp/2021029 www.mmnp-journal.org AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION Thomas Ashley1,*, Emilio Carrizosa1 and Enrique Fern´ andez-Cara2 Abstract. Dynamic optimisation provides complex challenges for optimal solution, but greatly increases applicability when considering time dependent situations. In this work, a constrained dynamic optimisation problem is analysed and subsequently applied to the resolution of a real-world engineering problem concerning Solar Power Tower plants. We study the existence of solutions and deduce an appropriate optimality characterisation in this applied framework. Two iterative algorithms are presented, convergence properties are discussed and a numerical illustration is given utilising realistic data. Finally, conclusions are drawn on the considered model and some ideas for future work are discussed. Mathematics Subject Classification. 90C46, 90C39, 37N40, 37M05. Received March 3, 2021. Accepted May 16, 2021. 1. Introduction Constrained optimisation problems are of great theoretical and practical interest, due to their complex nature and ability to describe real-world situations. Recent literature on the application of constrained optimisation can be seen from [13], where the authors consider on-line optimisation applied to a chemical reactor and [11], where a review on applied optimisation methods can be found. The mathematical complexity of such problems is amplified when the decision variables are functions of time and consequently belong to an infinite dimensional vector space. Related problems model many realistic situations in a very acceptable and efficient way. For example in [18], the authors utilise approximate dynamic programming to optimise problems within ambulance management and present a numerical illustration utilising real-world data; in [16], the authors present a review of dynamic optimisation vehicle routing problems and discuss various solution methods. In a previous work, the authors considered a stationary constrained problem motivated by the optimisation of the aiming action in a Solar Power Tower (SPT) plant, see [1]. The purpose of that paper was to determine strategies that maximise the radiation energy reaching the target (receiver). Optimisation techniques have been applied regularly in the literature to renewable energy research; see [8], where the authors present a review of related methods in sustainable buildings; see also [6], where multi-objective techniques are applied to Keywords and phrases: Dynamic optimisation, continuous optimisation, operations research, applications, aiming strategies for Solar Power Tower plants. 1IMUS – Instituto de Matem´aticas de la Universidad de Sevilla, Seville, Spain. 2Dep. EDAN and IMUS, Universidad de Sevilla, Seville, Spain. *Corresponding author: [email protected] c The authors. Published by EDP Sciences, 2021 This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 2T. ASHLEY ET AL. renewable energy systems, utilising particle swarm methods. The work presented in this paper looks to extend the optimisation model considered in [1] to the dynamic case, where additional time-dependent constraints are considered. In order to present a rigorous formulation of the problem we are interested in, let us introduce some notation and let us recall some definitions. As usual, for any Euclidean space S, we denote by |·|the associated norm. On the other hand, H1(0, T;S) stands for the Sobolev space of continuous functions p: [0, T ]7→ Ssuch that ˙ p(t) exists a.e. and ZT 0 |˙ p(t)|2dt < +∞. The associated duality product will be denoted by h· ,·i. Also, the symbol Cwill denote a generic positive constant. In this work, we consider extremal problems of the following kind: (Maximise J(p) = RT 0G(t, p(t)) dt Subject to p∈Pad,(1.1) where Pad is a subset of a Hilbert space Pof functions p=p(t) that take values in an Euclidean space F. It will be assumed that        G: [0, T]×F7→ Ris a continuous function, Gis differentiable with respect to pat any (t, p) and ∂G ∂p: [0, T]×F7→ Ris continuous. (1.2) More precisely, Pand Pad will be given by P=H1(0, T;F),(1.3) Pad ={p∈P:p(t)∈R∀t∈[0, T], M(p)≤σ},(1.4) where Ris a nonempty compact convex set in F, (1.5) M:H1(0, T;F)7→ Eis a C1mapping (Eis another Euclidean space) (1.6) and σ∈E. In (1.4) and henceforth, the inequality M(p)≤σmust be understood component-wise. In this paper, our aim is to solve theoretically and numerically the previous optimisation problem (1.1), paying special attention to applications to solar energy production. Details on the meaning of the variables and the functionals in the context of SPT plants will be given in Sections 4and 5. The plan of the paper is the following. In Section 2, we first prove the existence of a solution to (1.1) and subsequently characterise the solutions by an appropriate optimality system. We then formulate and analyse two iterative algorithms in Section 3, the first one relying on a penalty method and the second one using Augmented Lagrangian techniques. Section 4is devoted to particularise (1.1) in the context of a SPT plant; there, the existence and characterisation of optimal aiming strategies are established. The algorithms are illustrated with AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 3 numerical experiments for a real SPT plant in Section 5. Finally, Section 6contains some conclusions and the description of future work. 2. Theoretical analysis: existence and optimality results This section deals with the theoretical analysis of (1.1). To this purpose, let us introduce the set P0={p∈ P:p(t)∈R∀t∈[0, T]}. It will be assumed that the mapping Min (1.6) satisfies the following: Mis sequentially weakly lower semicontinuous in H1(0, T;F), (2.1) in the sense that pn→pweakly in H1(0, T;F) implies lim infn→+∞M(pn)≥M(p) and Mis bounded on bounded sets and coercive in P0, (2.2) in the sense that, for any σ∈E, the set of functions p∈P0satisfying M(p)≤σis bounded in H1(0, T;F) (equivalently, if pn∈P0for all nand kpnkp→+∞, then |(M(pn)−σ)+| → +∞for all σ∈E; here, z+denotes the positive part of z, understood component-wise). Our first result is the following: Theorem 2.1. Assume that the assumptions (1.2)–(2.2)are satisfied and the set Pad, given by (1.4), is nonempty. Then, there exists at least one solution to (1.1). Proof. The proof is standard if we take into account carefully the properties of the spaces and mappings involved in the formulation of the problem, see for instance [7]. Thus, let {pn}be a maximising sequence for (1.1), that is, a sequence in Pad such that J(pn)→sup p∈Pad J(p) as n→+∞. The pnare uniformly bounded in H1(0, T;F), since they all belong to Pad and (2.2) holds. Consequently, at least for a subsequence (again indexed by n), one has pn→ˆ pweakly in H1(0, T;F) and strongly in C0([0, T]; F). The set Pad is closed in P, thanks to (1.4), (1.5) and (2.1). Hence, ˆ p∈Pad. On the other hand, p7→ J(p) is continuous in C0([0, T ]; F), thanks to (1.2). Therefore, one has J(ˆ p) = max p∈Pad J(p). This ends the proof. Remark 2.2. In view of the possible non-convexity of Jand M, uniqueness is in general out of scope. From now on, it will be assumed that the hypotheses (1.2)–(1.6) are satisfied. Note that P0is a non-empty closed convex set of H1(0, T ;F). In our second result, we present suitable necessary optimality conditions that must be satisfied by the solutions to (1.1). For simplicity, we will assume from now on that E=R2and we will denote by M1and M2(resp. σ1and σ2) the components of M(resp. σ). Theorem 2.3. Let ˆ pbe a solution to (1.1). Assume that the constraints associated to Mare qualified at ˆp, that is: • ∃q1,q2∈H1(0, T;F)such that h∇Mi(ˆ p),qii<0for i= 1,2 4T. ASHLEY ET AL. Then, there exist λ1, λ2≥0such that the triplet (ˆ p, λ1, λ2)satisfies (h∇J(ˆ p),p−ˆ pi − λ1h∇M1(ˆ p),p−ˆ pi − λ2h∇M2(ˆ p),p−ˆ pi ≤ 0∀p∈P0, λi(Mi(ˆ p)−σi)=0, i = 1,2.(2.3) For the proof, we can use several arguments. Thus, we can directly apply the Karush-Kuhn-Tucker principle; see for instance Theorem 9-2.3 in [15]. We can also argue as in the proof of the Dubovirski-Milyouyin formalism; see [10] for details. 3. Some iterative algorithms 3.1. Penalisation In this section, we introduce an iterative method for the solution of (1.1) based on partial penalisation techniques. The advantage of this approach is that it reduces the task to the solution of other optimisation problems where only some constraints are kept (those easier to handle). The drawback is that a (small) parameter must be introduced and this can have a significant (undesired) influence at the numerical level. Thus, let us set Jµ(p) := J(p)−1 2µ|(M(p)−σ)+|2(3.1) for any µ > 0. Recall that (M(p)−σ)+stands for the positive part of M(p)−σand |·|denotes the Euclidean norm in E. Then, a suitable approximation of (1.1) is the following: (Maximise Jµ(p) Subject to p∈P0 (3.2) and it is reasonable to expect that, for small µ > 0, any solution to the penalized problem (3.2) solves approximately (1.1). Algorithm (Penalisation): 1. Choose µ0>0 and p0∈P0. 2. For n≥0, solve the following extremal problem:    Maximise ˆ Jn(p) := J(p)−1 2µn|(M(p)−σ)+|2 Subject to p∈P0 (3.3) 3. Take µn+1 =aµnwith 0 <a<1 and go to step 2 with nreplaced by n+ 1. Theorem 3.1. Let pnbe the solution to (3.3)for all n≥1, and assume that µn→0. Then, any weak limit p∗ of {pn}is a solution to (1.1). Proof. Let pbe a solution to (1.1). The following holds: 1. For all n≥0 there exists pn. Indeed, we are maximising in (3.3) a sequentially weakly upper semicontinuous function in a closed convex set P0⊂H1(0, T;F). Arguing as in the proof of Theorem 2.1, the existence of a solution follows. AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 5 2. The pnare uniformly bounded in H1(0, T;F), since J(pn)−1 2µn|(M(pn)−σ)+|2≥J(p)∀n≥1. Moreover, |(M(pn)−σ)+|2≤Cµnand limn→+∞|(M(pn)−σ)+|2= 0. So, there exists at least one weak limit point p∗. 3. If p∗is a weak limit, at least for a subsequence one has pk→p∗weakly in H1(0, T;F) and strongly in C0([0, T]; F). Then, |(M(p∗)−σ)+|2≤lim inf k→+∞|(M(pk)−σ)+|2= 0, whence pvec∗∈Pad. On the other hand, the following holds for any p0∈Pad: J(p0) = ˆ Jk(p0)≤ˆ Jk(pk) for all kwhence J(p0)≤lim inf k→+∞Jµk(pk)≤J(p∗). This ends the proof. From the practical viewpoint, instead of solving each subproblem (3.3) exactly, it can be more interesting to stop the resolution when the size of ∇ˆ Jnis sufficiently small. This motivates the following result: Theorem 3.2. Let us assume that µn&0,τn&0, |(M(p)−σ)+|2≥Cτnkpkp∀p∈P0 and lim inf n→+∞h∇M(pn),p0−pni ≤ h∇M(p∗),p0−p∗i ∀p0∈P0 and let us stop the resolution of (3.3)when k∇ ˆ Jn(p)kp=k∇J(p)−1 µn(M(p)−σ)+· ∇M(p)kp≤τn. Then, any weak limit p∗solves (1.1). Furthermore, if 1 µn(Mi(pn)−σi)+→λi ∗, i = 1,2 for a subsequence, then the Karush-Khun-Tucker conditions (2.3)are satisfied by p∗,λ1 ∗and λ2 ∗. Proof. Again, let pbe a solution to (1.1). Then: 1. If pndenotes the n-th iterate, with k∇ ˆ Jn(pn)kP≤τn, one has J(pn)−1 2µn|(M(pn)−σ)+|2≥J(p)−τnkpn−pkP. As a consequence, |(M(pn)−σ)+|2≤Cµn|(M(pn)−σ)+|2+ 1 6T. ASHLEY ET AL. and kpnkP≤C. We deduce that, at least for a subsequence, pk→p∗weakly in H1(0, T;F) and strongly in C0([0, T]; F). 2. Any weak limit p∗is a solution to (1.1) (proved as in Thm. 3.1). Now, let us set λi n=1 µn(Mi(pn)−σi)+ for i= 1,2. Then, |λ1 n|2+|λ2 n|2=2 µn(1 µn|(M(pn)−σ)+|2)≤C[|(M(pn)−σ)+|2+ 1] ≤C. Therefore, we can suppose that λi n→λi ∗, where λi ∗≥0. Since k∇ ˆ Jn(pn)kP=k∇J(pn)−1 µn 2 X i=1 (Mi(pn)−σi)∇Mi(pn)kP =k∇J(pn)− 2 X i=1 λi n∇Mi(pn)k ≤ τn, one has h∇J(pn),pni → h∇J(p∗),p∗iand h∇J(p∗),p0−pni → h∇J(p∗),p0−p∗i ∀p0∈P. 3. Consequently, for any p0∈P0, we see that h∇J(p∗)− 2 X i=1 λi ∗∇Mi(p∗),p0−pni =h∇J(pn)− 2 X i=1 λi n∇Mi(pn),p0−pni+h∇J(p∗)− ∇J(pn),p0−pni + 2 X i=1 (λi n−λi ∗)h∇Mi(pn),p0−pni+ 2 X i=1 λi ∗h∇Mi(pn)− ∇Mi(p∗),p0−pni. This yields h∇J(p∗)− 2 X i=1 λi ∗∇Mi(p∗),p0−pni ≤τnkp0−pnkP+|h∇J(p∗)− ∇J(pn),p0−pni| + 2 X i=1 |h∇Mi(pn),p0−pni|λi n−λi ∗| + 2 X i=1 λi ∗(h∇Mi(pn),p0−pni− <∇Mi(p∗),p0−pn>) and, taking lower limits, we easily find that h∇J(p∗)−P2 i=1 λi ∗∇Mi(p∗),p0−p∗i ≤ 0 for all p0∈P0. This ends the proof. In the remainder of this section, we will be concerned with the numerical solution of (3.2). To this purpose, we will first introduce a time discretisation and then an iterative algorithm of the gradient kind. Thus, let us begin by replacing H1(0, T ;F) by a finite dimensional space. The easiest and most natural way is to introduce a large integer N, set τ:= T/N, consider a uniform partition {t0= 0 < t1< t2< ... < tN=T} AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 7 with tn=nτ for all nand, then, work in the corresponding space of continuous and piecewise linear functions pN: [0, T]7→ F. This space will be denoted by XN. We will also consider the closed convex set P0,N :=XN∩P0and the orthogonal projector P0,N :XN7→ P0,N . Observe that P0,N (pN)(t) is, for each t, the projection of each component of pN(t) onto R. For fixed µ > 0, the considered N-th approximated problem to (3.2) is as follows: (Maximise Jµ(pN) Subject to pN∈P0,N .(3.4) In order to solve this problem, we can apply a gradient ascent algorithm with variable step size and projection. Accordingly, in the n-th iterate we compute the function pn+1 N, with pn+1 N=P0,N (˜ pn+1 N),˜ pn+1 N=pn N+γn∇Jµ(pn N).(3.5) Here, γnis a conveniently chosen positive number and ∇Jµ(pn N) denotes the gradient of the objective function in (3.4), that is, ∇Jµ(pN) = ∇J(pN)−1 µ(M(pN)−σ)+· ∇M(pN).(3.6) Therefore, the gradient ascent algorithm requires, at each step, the calculation of the derivative of the objective function and a projection of each component of pN(t) for each nodal time t=tj. Obviously, the complexity of this computation depends on dim F,Nand the properties of the particular function G=G(t, p) and the mapping M:H1(0, T;F)7→ E. The previous penalisation algorithm, denoted Algorithm 1in this paper, is outlined in a pseudocode shown below. Algorithm 1 Penalisation µ←Set penalisation parameter (small) float γ0←Set initial step size float Tol ←Simulation tolerance float p0←Initial variable set float for N= 1 to T do while ObjDiff >Tol do while StepDiff ≤0do J(p)←Calculate objective function float M(p)←Calculate constraint set float ∇Jµ(pN) = ∇J(pN)−1 µ(M(pN)−σ)+· ∇M(pN). ˜ pn+1 N=pn N+γn∇Jµ(pn N). pn+1 N=P0,N (˜ pn+1 N) J(˜ pn+1 N)←Calculate updated objective function float StepDiff =|Jµ(pn+1 N)−Jµ(pn N)| γn+1 =γn· Next; ObjDiff =|Jµ(pn+1 N)−Jµ(pn+1 N−1)| Next; end; 8T. ASHLEY ET AL. 3.2. Augmented Lagrangian The optimisation problem (1.1) can also be solved using the information furnished by Theorem 2.3. This is performed in this section using Augmented Lagrangian techniques. Again, it will be assumed that E=R2and M1and M2will denote the components of M. Thus, let us first introduce the so called Augmented Lagrangian Lµ:H1(0, T;F)×E7→ R, with Lµ(p;λ) := J(p)− 2 X i=1 ψ(Mi(p)−σi, λi;µ),(3.7) where we have set ψ(z, β;µ) :=    z·β+1 2µ|z|2if z+µβ ≥0 −µ 2|β|2otherwise (3.8) and, again, µ > 0. Then, it can be proved that the (1.1) is equivalent to the following extremal problem    Minimise sup p∈P0 Lµ(p;λ) Subject to: λ∈E, λ ≥0. (3.9) Let us provide an explanation. The inequality constraint M(p)≤σis equivalent to the equality M(p)+s=σ, with the so called slack variable sbelonging to Eand satisfying s≥0. In view of the optimality conditions (2.3), it makes sense to consider the “Modified Lagrangian” ˜ L(p, s;λ) := J(p)−λ·(M(p) + s−σ) (3.10) and its penalised version ˜ Lµ(p, s;λ) := J(p)−λ·(M(p) + s−σ)−1 2µ|M(p) + s−σ|2.(3.11) Then, recalling (3.7)–(3.8), it is easy to check that sup s∈E, s≥0 ˜ Lµ(p, s;λ) = Lµ(p;λ), whence we see that (3.9) is an appropriate reformulation of (1.1). Algorithm (Augmented Lagrangian): 1. Choose µ0, λ0 1, λ0 2≥0 adn p0∈P0. 2. For n≥0, solve the following extremal problem:      Maximise ˆ Ln(p;λn) := J(p)− 2 X i=1 ψ(Mi(p)−σi, λn i;µn) Subject to p∈P0. (3.12) AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 9 3. Update the multipliers: λn+1 i= (λn i+1 µn(Mi(pn)−σi))+. 4. Take µn+1 =aµnwith 0 <a<1 and go to step 2 with nreplaced by n+ 1. The following results justify the use of the augmented Lagrangian Lµin (3.7) for small values of µ: Theorem 3.3. Assume that Jand the constraint functionals Miare in C2in Pand p∗is a solution to (1.1) (or at least a local maximiser), with M(p∗)< σ. Also, assume that the usual second order sufficient conditions hold at p∗for some λ∗= (λ1 ∗, λ2 ∗)with λi ∗≥0. Then, there exists ¯µ > 0such that, for all µ∈(0,¯µ],p∗is a minimiser of p→ Lµ(p;λ∗). Proof. First, note that for small ¯µand 0 <µ<¯µ,ψ(Mi(p∗)−σi, λi ∗;µ)≡ −1 2µ|λi ∗|2, whence ∇pLµ(p;λ∗)|p=p∗=∇J(p∗), and L00 µ(p;λ∗)|p=p∗=J00(p∗). By assumption, p∗solves (1.1), Mi(p∗)> σi, h∇J(p∗)− 2 X i=1 λi ∗∇Mi(p∗),p0−p∗i ≤ 0∀p0∈P0, λi ∗≥0 and λi ∗(Mi(p∗)−σi) = 0 (whence λi ∗= 0) and (J00 − 2 X i=1 λi ∗M00 i)(p∗;p0−p∗,p0−p∗)≤ −akp0−p∗k2 P∀p0∈P0(a > 0). Consequently, the second order sufficient conditions hold for p→ Lµ(p;λ∗) at p∗if µ≤¯µand the proof is done. Theorem 3.4. Assume that the hypotheses in Theorem 3.3 are satisfied with λ1 ∗=λ2 ∗= 0. Then, there exists a, ¯µ,  > 0such that, if |λ1|+|λ2| ≤ a/µ and 0< µ ≤¯µ, one has: 1. The extremal problem      Maximise Lµ(p;λ) := J(p)−P2 i=1 ψ(Mi(p)−σi, λi;µ) Subject to p∈P0,kp−p∗kP≤ kp−p∗kP≤Cµ(|λ1|+|λ2|) (3.120) is uniquely solvable. 2. Let ˜ pbe the unique solution to (3.120). Then hL00 µ(˜ p;λ)(p0−˜ p),p0−˜ pi≤−˜akp0−˜ pk2 P∀p0∈P0.(3.13) Proof. From standard arguments, we readily see that (3.120) is solvable. On the other hand, for small ¯µand , one has ∇Lµ(p;λ) = ∇J(p) and hL00 µ(p;λ)q, zi=J00(p;q, z)∀q, z ∈P. Consequently, one has (3.13) at any solution ˜ p. This implies uniqueness and the proof is done. 16 T. ASHLEY ET AL. It is not difficult to compute from this identity the partial derivatives of JN µand, accordingly, its gradient ∇Jµ. This makes it possible to apply the gradient ascent algorithm in the present context. The choice of step size parameters γkused in the iterates (3.5) have a large impact on the speed of convergence of the algorithm. Much research has been conducted in the optimisation of this parameter, see for example [12]. A frequently successful strategy corresponds to the so called Armijo’s Rule [14]. In this method, a constant value ∈[0,1] is used to iteratively reduce γkuntil an increase in the objective function is no longer found. Thus, with a suitable choice of an initial γ0, the largest possible value is taken per iteration. In this work, we apply a variant to the Armijo’s Rule, whereby the parameter γkis different for each component, that is, for each heliostat h∈H. Consequently, the γk hare computed according to the rules γk h=γk−1 h·h, h = 1, . . . , H. (4.15) The gradient ascent method with projection is then: pn+1 N=P0,N (˜ pn+1 N),˜ pn+1 N=pn N+ Γk∇Jµ(pn N),(4.16) where Γk= diag(γk 1, . . . , γk N). The gradient ascent method applied to the penalisation algorithm, as described in the pseudocode given in Algorithm 2, is applied at each step of the discretised time period. The computed optimal result at a given time is then used as the initial guess for the next time step. In this way, the overall optimal schedule can be found by considering each discretised time point individually. 4.2.2. Augmented Lagrangian algorithm Following the finite dimensional reduction method outlined in Section 3, the Augmented Lagrangian given in equations (3.7) can be approximated as follows: Lµ(pN, λ)≃ LN µ(pN, λ):=A N X n=1 I X i=1 H X h=1 Fui,vi(h, ph(tn), tn) −(1−A) N X n=1 I X i=1 H X h=1 b(tn,p(tn); ui, vi)2− 2 X i=1 ψ(MN i(pN)−σi, λi;µ), (4.17) where ψis defined in (3.8) and MN 1and MN 2are respectively given in (4.13) and (4.14). The gradient of this function can be calculated as before and the gradient ascent method can be applied to the intermediate problems (3.15). The Augmented Lagrangian algorithm, as summarised in the pseudocode in Section 3.2, follows a similar procedure to that of the penalisation algorithm, with extra stages to update the Lagrangian multipliers. Recall that, after the computation of the solution to (3.15), the Lagrange multipliers must be updated, according to the formula (3.150). Then, a new optimisation problem must be solved using the last computed solution as an initial iterate and so on. In the following section, we present an illustrative example that demonstrates, in the framework of the optimisation of SPT plants, the usefulness of the formulation (1.1) and the theoretical results in Section 2and the functionality of the previous algorithms. 5. A numerical experiment The behavior of the previous algorithms is illustrated for the PS10 SPT plant located in Sanl´ucar la Mayor, Seville (Spain); see http://www.abengoa.com/web/en. The field of this SPT plant has 624 heliostats arranged AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 17 Figure 1. Heliostat layout. Figure 2. Midday radiation. in a radial pattern around a centrally located tower. The layout of the heliostats can be seen in Figure 1, where the receiver is mounted atop a North facing tower. In Figure 2, the reflected solar radiation at midday is shown, where the heliostats are colour coded according to the energy they would provide to the system if they aimed at the centre of the receiver. From this figure, it can be clearly seen that an adequate aiming strategy is important, as there are large differences in the energy contributions provided by the heliostats in the field. It is completely natural to fix a dynamic aiming strategy, as the incident radiation on the heliostat field changes over time. Utilising the algorithms developed in Section 3and the SPT plant model described in Section 4, we present a numerical illustration that finds the optimal aiming strategy for the PS10 SPT plant across a single day. We consider 10 equally spaced time points, where input to the model includes incident solar radiation on the field and solar angle. We look to optimise the general dynamic optimisation problem (1.1), by considering the two forms given by the algorithms: penalisation (3.4) and Augmented Lagrangian (3.14). As already mentioned, the objective 18 T. ASHLEY ET AL. Table 1. Penalisation algorithm parameter values. Parameter Value Summary T10 Number of time points σ10M1constraint threshold σ21e4 M2constraint threshold H624 Number of heliostats in field. µ1e5 Penalisation constant Vp0.5 Velocity constraint constant Etar 2.2e6 Target uniform distribution value (Wm2) A0.7 Weighting parameter γ00.01 Initial step size is to maximise radiation reaching the receiver surface across the day, whilst restricting the movement speed of heliostat aim points and also limiting the change in radiation over time at any point on the receiver. 5.1. Penalisation algorithm The penalisation algorithm is applied to the 10 point optimisation problem, with the parameter values given in Table 1. We set a uniform target distribution and fix limits on the aiming point velocities and global change in radiation along the receiver. Considering the constraints derived in Section 3.1, we then look to maximise the objective function at each time point utilising the gradient ascent method. To start the algorithm, we define an initial set of aiming points on the receiver, randomly spread, as shown in Figure 3. This choice is ideal for early morning, as it allows for a slow warm-up of the entire receiver surface. The resultant set of aiming points and the radiation distribution on the receiver are then used as the initial solution to the first time problem. Each subsequent time point is then considered, utilising the computed optimal solution from the previous one as its initial solution. The aiming strategies for all time points are given in Figures 4–13. From Figures 4–13, the evolution of the aiming strategy across the day can be seen. The effect of the dynamic constraints on aiming point velocity and change in radiation can be seen in the slow movement of aim points towards the center in Figures 4–8. During these early hours of low incident solar radiation, the maximum is captured when a heliostat aims at the center. The change in aiming strategy from Figures 8and 9identifies a shift in behaviour, caused by the incident radiation on the heliostat field and the target distribution imposed in the constraints. Figure 14 gives the level of incident radiation on the heliostat field for the considered time period, and Figure 15 gives the maximum radiation value detected on the receiver with the indicated aiming strategy. It can be seen that, as the level of incident radiation on the field increases (towards midday), the same happens to the maximum level of radiation on a particular point on the receiver. This is due to the centrally focused aiming strategy shown in Figure 7. Once the level of radiation reaches the target distribution limit given in Table 1, at time t= 5, the aiming strategy must adjust in order to maintain a uniform distribution, as seen in Figure 9. The computation time of the penalisation algorithm is dependent on the prescribed model parameters, such as step size, constraint limits, convergence tolerance, etc. However, with an adequate selection of these parameters, the previous numerical illustration can be achieved in 2.5 minutes utilising a computer with specifications: Intel R CoreTM i7-7700HQ CPU @ 2.80 GHz. 5.2. Augmented Lagrangian algorithm The Augmented Lagrangian algorithm is also applied to the 10 point optimisation problem, with the parameter values given in Table 2. We set a uniform target distribution and fix limits on the aiming point velocities AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 19 Figure 3. Initial aiming point distribution. Figure 4. t= 1 aiming distribution. 20 T. ASHLEY ET AL. Figure 5. t= 2 aiming distribution. Figure 6. t= 3 aiming distribution. AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 21 Figure 7. t= 4 aiming distribution. Figure 8. t= 5 aiming distribution. 22 T. ASHLEY ET AL. Figure 9. t= 6 aiming distribution. Figure 10. t= 7 aiming distribution. AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 23 Figure 11. t= 8 aiming distribution. Figure 12. t= 9 aiming distribution. 24 T. ASHLEY ET AL. Figure 13. t= 10 aiming distribution. Figure 14. Incident radiation on SPT plant. and change in radiation along the receiver. Considering the constraints in Section 3.2, we then look to maximise the objective function across each time point utilising the gradient ascent method. To start the algorithm, we define an initial set of aiming points on the receiver as in the case of the penalisation algorithm, shown in Figure 16. For each time point considered, the algorithm takes the initial solution from the previous time step, and optimises using the gradient ascent method described in Section 3.2. Once a solution has been found for a AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 25 Figure 15. Maximum flux on receiver. Table 2. Augmented Lagrangian algorithm parameter values. Parameter Value Summary T10 Number of time points σ10M1constraint threshold σ21e4 M2constraint threshold H624 Number of heliostats in field. µ2e−5 Penalisation constant λ01e6 Initial Lagrange multiplier Vp0.5 Velocity constraint constant Etar 2.2e6 Target uniform distribution value (Wm2) A0.7 Weighting parameter γ00.01 Initial step size particular time point, the Lagrangian multipliers are updated, and the solution is used to re-optimise the same time point. This is repeated until no improvement on the solution can be found. The resultant solution is then used as the initial solution to the next time point. The computed aiming strategies using the augmented lagrangian algorithm are then given in Figures 17–41. The maximum flux on the receiver surface and the incident radiation level on the heliostat field are given in Figure 42. The results furnished by the Augmented Lagrangian algorithm agree with those provided by the penalisation algorithm. 5.3. Numerical considerations The PS10 plant field contains 624 heliostats, however there exist other larger SPT plants (with potentially thousands of heliostats) which greatly increases the dimensionality of the problem. The aiming points are described on the receiver surface in polar coordinates, necessitating 2Hvariables for the gradient ascent technique. Within each step of the algorithm, the gradients must be calculated, and aiming points updated, before the radiation calculations are performed. For larger values of Hincreases to sufficient levels, the application of the algorithms presented in this paper may become infeasible, due to computation times. For the numerical approximations presented in this paper, the heliostat field size will have an effect on performance. 32 T. ASHLEY ET AL. Figure 28. t= 4, λ1aiming distribution. Figure 29. t= 5, λ0aiming distribution. AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 33 Figure 30. t= 5, λ1aiming distribution. Figure 31. t= 5, λ2aiming distribution. 34 T. ASHLEY ET AL. Figure 32. t= 6, λ0aiming distribution. Figure 33. t= 6, λ1aiming distribution. AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 35 Figure 34. t= 7, λ0aiming distribution. Figure 35. t= 7, λ1aiming distribution. 36 T. ASHLEY ET AL. Figure 36. t= 8, λ0aiming distribution. Figure 37. t= 8, λ1aiming distribution. AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 37 Figure 38. t= 9, λ0aiming distribution. Figure 39. t= 9, λ1aiming distribution. 38 T. ASHLEY ET AL. Figure 40. t= 10, λ0aiming distribution. Figure 41. t= 10, λ1aiming distribution. AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION 39 Figure 42. Maximum flux on receiver and incident solar radiation. 6. Conclusions and future work In this work, a general dynamic optimisation problem has been investigated. Theoretical properties of the problem have been discussed, including the demonstration of solution existence and optimality. Furthermore, a real-world example for a SPT plant has been considered, where a mathematical model that describes transfer of energy within the system is developed. The theoretical properties of the general optimisation problem are shown applicable to the real world problem, and dynamic constraints that describe SPT plant limitations are given. Two algorithms have been considered to find the optimal solution, and a numerical illustration is given using real data. The numerical illustration finds the optimal aiming strategy for a SPT plant over the period of one day, considering the change in incident radiation from the Sun as input. The physical limitations of the plant have been introduced as dynamic constraints, in terms of heliostat rotation speed and flux homogeneity on the receiver surface over time. The two algorithms presented in this paper provide similar results in the numerical experiment, in similar computational times. The Augmented Lagrangian algorithm is a modification of a penalty technique, which should provide better numerical stability in some cases. With larger size problems, and adequately chosen parameters, this model could increase the reliability of the algorithm. The dynamic optimisation problem considered in this paper must be viewed as an improvement of other research concerning the optimisation of aiming strategies in SPT plants, for example [4,19], due to the inclusion of dynamic constraints. Instead of optimising the aiming strategy at certain fixed times, the method presented in this paper looks to optimise across a time period. This approach can arrive to a solution that more closely reflects the true optimum when considering problems in dynamic systems. The methods presented in this work can be adapted to all types of SPT plants, and even other forms of concentrating solar power technology. The inclusion of more heliostats, or multiple receivers, has been carried out in several real-world plants and the work presented here can be directly extended to consider these cases. The use of a modified Armijos’ Rule for the step size in the gradient ascent method can lead to faster convergence. However, with a poor parameter selection, this may actually reduce the convergence speed. Further adaptation to the algorithms presented in this work should consider carefully the effects of such techniques, in conjunction with a highly multi-modal objective function. Considering the problem dimensionality when increasing the number of variables, specifically the number of heliostats in our example, one method that could be used to extend the presented approach could be to use a 40 T. ASHLEY ET AL. clustering algorithm. Thus, the heliostats could be clustered into groups, considering a weighted objective of difference in location and energy, as shown in [3]. An extension to this work would be the integration of a thermal transfer model of the receiver, whereby the aiming strategy and thermal transfer could be coupled within the optimisation process. This would further increase the level to which the true system is modelled, as dynamic transfer of heat through the receiver changes based upon flow rates within the tubes. This would therefore have an affect on the aiming strategy. Acknowledgements. This research is being supported by the Spanish Government through the research project PCIN- 2015-108, and in part through the EC H2020 MSCA RISE NeEDS Project (Grant agreement ID: 822214) and is being conducted within the group MTM2015-65915-R at the University of Seville. References [1] T. Ashley, E. Carrizosa and E. Fern´andez-Cara, Optimisation of aiming strategies in Solar Power Tower plants. Energy (2017). [2] T. Ashley, E. Carrizosa and E. Fern´andez-Cara, Continuous Optimisation Techniques for Optimal Aiming Strategies in Solar Power Tower Plants (Submitted) (2018). [3] T. Ashley, E. Carrizosa and E. Fern´andez-Cara, Heliostat field cleaning scheduling for Solar Power Tower plants: a heuristic approach. Appl. Energy 235 (2019) 653–660. [4] M. Astolfi, M. Binotti, S. Mazzola, L. Zanellato and G. Manzolini, Heliostat aiming point optimization for external tower receiver. Solar Energy 157 (2017) 1114–1129. [5] R. Braun, C. Diep, W.T. Hamilton, M.J. Wagner, J. Dent and A. Newman, Optimizing dispatch for a concentrated solar power tower. Solar Energy 174 (2018) 1198–1211. [6] D.P. Clarke, Y.M. Al-Abdeli and G. Kothapalli, Multi-objective optimisation of renewable hybrid energy systems with desalination. Energy 88 (2015) 457–468. [7] I. Ekeland and R. Temam, Convex analysis and variational problems. Elsevier (1976). [8] R. Evins, A review of computational optimisation methods applied to sustainable building design. Renew. Sustain. Energy Rev. 22 (2013) 230–245. [9] O. Farges, J.J. B´ezian and M. El Hafi, Global optimization of solar power tower systems using a Monte Carlo algorithm: application to a redesign of the PS10 solar thermal power plant. Renew. Energy 119 (2018) 345–353. [10] I.V. Girsanov, Lectures on mathematical theory of extremum problems. Springer-Verlag, Berlin-New York (1972). [11] F. Kheiri, A review on optimization methods applied in energy-efficient building geometry and envelope design. Renew. Sustain. Energy Rev. 92 (2018) 897–920. [12] Z. Liu and H. Liu, Several efficient gradient methods with approximate optimal stepsizes for large scale unconstrained optimization. J. Comput. Appl. Math. 328 (2018) 400–413. [13] M. Mansour and J.E. Ellis, Methodology of on-line optimization applied to a chemical reactor. Appl. Math. Model. 32 (2008) 170–184. [14] G.P. McCormick, A modification of Armijo’s step-size rule for negative curvature. Math. Program. 13 (1977) 111–115. [15] J. Nocedal, S.J. Wright and S.M. Robinson, Numerical Optimization, Springer, 1 ed. (1999). [16] V. Pillac, M. Gendreau, C. Gu´eret and A.L. Medaglia, A review of dynamic vehicle routing problems. Eur. J. Operat. Res. 225 (2013) 1–11. [17] S. Relloso and E. Garc´ıa, Tower technology cost reduction approach after Gemasolar experience. Energy Proc. 69 (2015) 1660–1666. [18] V. Schmid, Solving the dynamic ambulance relocation and dispatching problem using approximate dynamic programming. Eur. J. Oper. Res. 219 (2012) 611–621. [19] K. Wang, Y.-l. He, X.-d. Xue and B.-c. Du, Multi-objective optimization of the aiming strategy for the solar power tower with a cavity receiver by using the non-dominated sorting genetic algorithm. Appl. Energy 205 (2017) 399–416. [20] Q. Yu, Z. Wang and E. Xu, Analysis and improvement of solar flux distribution inside a cavity receiver based on multi-focal points of heliostat field. Appl. Energy 136 (2014) 417–430.