AN EXTREMAL PROBLEM WITH APPLICATIONS TO RENEWABLE ENERGY PRODUCTION

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.


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.
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, H 1 (0, T ; S) stands for the Sobolev space of continuous functions p : [0, T ] → S such thaṫ p(t) exists a.e. and T 0 |ṗ(t)| 2 dt < +∞.
The associated duality product will be denoted by · , · . Also, the symbol C will denote a generic positive constant.
In this work, we consider extremal problems of the following kind: Maximise J(p) = T 0 G(t, p(t)) dt Subject to p ∈ P ad , (1.1) where P ad is a subset of a Hilbert space P of functions p = p(t) that take values in an Euclidean space F. It will be assumed that G is differentiable with respect to p at any (t, p) and ∂G ∂p : [0, T ] × F → R is continuous. (1.2) More precisely, P and P ad will be given by where R is a nonempty compact convex set in F, (1.5) M : H 1 (0, T ; F) → E is a C 1 mapping (E is 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 4 and 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 4 is 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 numerical experiments for a real SPT plant in Section 5. Finally, Section 6 contains some conclusions and the description of future work.

Theoretical analysis: existence and optimality results
This section deals with the theoretical analysis of (1.1). To this purpose, let us introduce the set P 0 = {p ∈ P : p(t) ∈ R ∀t ∈ [0, T ]}. It will be assumed that the mapping M in (1.6) satisfies the following: M is sequentially weakly lower semicontinuous in H 1 (0, T ; F), (2.1) in the sense that p n → p weakly in M is bounded on bounded sets and coercive in P 0 , (2.2) in the sense that, for any σ ∈ E, the set of functions p ∈ P 0 satisfying M (p) ≤ σ is bounded in H 1 (0, T ; F) (equivalently, if p n ∈ P 0 for all n and p n p → +∞, then |(M (p n ) − σ) + | → +∞ for all σ ∈ E; here, z + denotes the positive part of z, understood component-wise).
Our first result is the following: 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 {p n } be a maximising sequence for (1.1), that is, a sequence in P ad such that The p n are uniformly bounded in H 1 (0, T ; F), since they all belong to P ad and (2.2) holds. Consequently, at least for a subsequence (again indexed by n), one has p n →p weakly in H 1 (0, T ; F) and strongly in C 0 ([0, T ]; F).
The set P ad is closed in P, thanks to (1.4), (1.5) and (2.1). Hence,p ∈ P ad . On the other hand, p → J(p) is continuous in C 0 ([0, T ]; F), thanks to (1.2). Therefore, one has This ends the proof. From now on, it will be assumed that the hypotheses (1.2)-(1.6) are satisfied. Note that P 0 is a non-empty closed convex set of H 1 (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 = R 2 and we will denote by M 1 and M 2 (resp. σ 1 and σ 2 ) the components of M (resp. σ).
Theorem 2.3. Letp be a solution to (1.1). Assume that the constraints associated to M are qualified atp, that is: Then, there exist λ 1 , λ 2 ≥ 0 such that the triplet (p, λ 1 , λ 2 ) satisfies 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.

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 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: and it is reasonable to expect that, for small µ > 0, any solution to the penalized problem (3.2) solves approximately (1.1).
Theorem 3.1. Let p n be the solution to (3.3) for all n ≥ 1, and assume that µ n → 0. Then, any weak limit p * of {p n }is a solution to (1.1).
Proof. Let p be a solution to (1.1). The following holds: 1. For all n ≥ 0 there exists p n . Indeed, we are maximising in (3.3) a sequentially weakly upper semicontinuous function in a closed convex set P 0 ⊂ H 1 (0, T ; F). Arguing as in the proof of Theorem 2.1, the existence of a solution follows.
2. The p n are uniformly bounded in H 1 (0, T ; F), since Moreover, |(M (p n ) − σ) + | 2 ≤ Cµ n and lim n→+∞ |(M (p n ) − σ) + | 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 p k → p * weakly in H 1 (0, T ; F) and strongly in C 0 ([0, T ]; F). Then, whence pvec * ∈ P ad . On the other hand, the following holds for any p ∈ P ad : 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 ∇Ĵ n is sufficiently small. This motivates the following result: Then, any weak limit p * solves (1.1). Furthermore, if for a subsequence, then the Karush-Khun-Tucker conditions (2.3) are satisfied by p * , λ 1 * and λ 2 * . Proof. Again, let p be a solution to (1.1). Then: 1. If p n denotes the n-th iterate, with ∇Ĵ n (p n ) P ≤ τ n , one has As a consequence, |(M (p n ) − σ) + | 2 ≤ Cµ n |(M (p n ) − σ) + | 2 + 1 and p n P ≤ C. We deduce that, at least for a subsequence, p k → p * weakly in H 1 (0, T ; F) and strongly in C 0 ([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 (M i (p n ) − σ i ) + for i = 1, 2. Then, Therefore, we can suppose that one has ∇J(p n ), p n → ∇J(p * ), p * and 3. Consequently, for any p ∈ P 0 , we see that

This yields
and, taking lower limits, we easily find that 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 H 1 (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 {t 0 = 0 < t 1 < t 2 < ... < t N = T } with t n = nτ for all n and, then, work in the corresponding space of continuous and piecewise linear functions This space will be denoted by X N . We will also consider the closed convex set P 0,N := X N ∩ P 0 and the orthogonal projector P 0,N : X N → P 0,N . Observe that P 0,N (p N )(t) is, for each t, the projection of each component of p N (t) onto R. For fixed µ > 0, the considered N -th approximated problem to (3.2) is as follows: 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 p n+1 N , with Here, γ n is a conveniently chosen positive number and ∇J µ (p n N ) denotes the gradient of the objective function in (3.4), that is, 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 p N (t) for each nodal time t = t j . Obviously, the complexity of this computation depends on dim F, N and the properties of the particular function G = G(t, p) and the mapping M : The previous penalisation algorithm, denoted Algorithm 1 in this paper, is outlined in a pseudocode shown below.

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 = R 2 and M 1 and M 2 will denote the components of M .
Thus, let us first introduce the so called Augmented Lagrangian L µ : where we have set and, again, µ > 0. Then, it can be proved that the (1.1) is equivalent to the following extremal problem Let us provide an explanation. The inequality constraint M (p) ≤ σ is equivalent to the equality M (p) + s = σ, with the so called slack variable s belonging to E and satisfying s ≥ 0. In view of the optimality conditions (2.3), it makes sense to consider the "Modified Lagrangian" and its penalised versionL Then, recalling (3.7)-(3.8), it is easy to check that sup s∈E, s≥0L whence we see that (3.9) is an appropriate reformulation of (1.1).
The following results justify the use of the augmented Lagrangian L µ in (3.7) for small values of µ: Theorem 3.3. Assume that J and the constraint functionals M i are in C 2 in P and 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 Proof. First, note that for smallμ and 0 < µ <μ, By assumption, Consequently, the second order sufficient conditions hold for p → L µ (p; λ * ) at p * if µ ≤μ and the proof is done.
Consequently, one has (3.13) at any solutionp. This implies uniqueness and the proof is done.
Arguing as in [15], it can also be proved that p − p * P ≤ Cµ(|λ 1 | + |λ 2 |) for some C > 0. As in the previous section, in order to solve (3.9) we must provide a finite dimensional approximation. With the notation used in Section 3.1, a suitable choice is the following: Subject to: λ ∈ E, λ ≥ 0. (3.14) This problem can be solved with a duality-penalty algorithm that at the nth step furnishes the multipliers λ n+1 i . The strategy is as follows: • Compute a solution p n N to the problem • Then, take As in Section 3.1, the solution of (3.15) can be carried out through a variable step gradient ascent method with projection. Thus, at the k-th step, we compute p n,k+1 N , with Algorithm 2 is outlined in the pseudocode shown below. In the following section, we will detail the model for a SPT plant and verify that all assumptions made in the formulation of the general problem in Section 1 remain valid.

The model for a SPT plant
The extremal problem (1.1) can be used to model the optimisation of the aiming strategy in a SPT plant if we assume that • p = p(t) defines the set of points aimed by the heliostats on the receiver at times t ∈ [0, T ], • R is the receiver surface, • G = A · G 1 + (1 − A) · G 2 is a balanced combination of the objective functions G 1 and G 2 , respectively corresponding to radiation and closeness to a target aiming strategy), • J(p) is, accordingly, a quantification of the payoff produced by the aiming strategy determined by p and • M (p) is a measure of the change over time of the aiming strategy, corresponding to p and the associated energy reaching the receiver.
The ith test point on the receiver surface A SPT plant is a renewable energy system, containing a field of heliostats with a centrally located tower, atop of which one has a receiver. The concentration of incident solar radiation onto the receiver by the heliostats allows for high temperatures to be achieved, which can then be used to drive a turbine through the use of a heat transfer medium.
Much research considering this form of renewable energy generation can be found in the literature, including [5], where the authors look to optimise annual revenue by considering dispatch measures. Also, the redesign of a known SPT plant heliostat field through the use of optimisation techniques is considered in [9].
As considered by the authors in previous work, see [1,2], the aiming point of each heliostat in the field on the receiver surface directly affects the energy generation. In particular, it has been shown that, although maximising incident radiation is important, it is also beneficial to maintain a desired flux distribution to aid thermal transfer [17,20].
In this paper we extend the models to the dynamic case, considering time dependent variables and constraints. These constraints, as introduced in Sections 1 and 2, concern physical aspects of the SPT. One of them must be viewed as a limitation on the rotational speed of the heliostat and is obviously justified by operational reasons. On the other hand, the radiation at any point of the receiver surface should not vary drastically over short time periods, in order to prevent flash heating and this is considered in the second constraint.
As in [1], the radiation passing through the system is modelled using a Gaussian distribution. We assume that there are H heliostats in the field and we denote by p h (t) the point aimed by the h-th heliostat at time t. We have p h (t) ∈Ω for all h and t; Ω ⊂ R 2 is bounded, open and convex with 0 ∈ Ω and we take R =Ω. Accordingly, the Euclidean space F will have dimension 2H: (4.1) The usual Euclidean norm in F will be denoted by · .
At time t, the radiation measured at a point (u, v) ∈ R and furnished by the h-th heliostat is given by F u,v (h, p h (t), t). For each h, it will be assumed that the real-valued function (u, v, p h , t) → F u,v (h, p h , t) is smooth. The associated total radiation captured over the receiver surface R for a given heliostat h is thus given by Therefore, the total radiation supplied by all heliostats at time t can be written in the form where we have used the notation p := (p 1 , . . . , p H ). In order to limit the motion of heliostats, so that they do not move faster than their velocity limits, we introduce the velocityṗ h =ṗ h (t) of each heliostat h, the velocity vectorṗ =ṗ(t) and a target velocity vector V p ∈ F with V p ≥ 0 and we impose the collective velocities to be approximately below V p in the sense that for some σ 1 ≥ 0. The change in time in radiation at each point of the receiver must also be limited and kept below a fixed value σ 2 . This will be witten in the form Our objective is to maximise the radiation captured by the receiver, whilst maintaining a target energy profile which maximises absorption. Thus, we deal with the constrained optimisation problem (1.1), where the σ i are prescribed non-negative constants and E u,v tar = E u,v tar (t) is the desired target distribution at time t. For convenience, we will assume that This is sufficient to guarantee that the zero function in H 1 (0, T ; F) belongs to P ad and, consequently, P ad is non-empty. From a realistic viewpoint, the assumption on the smoothness of (u, v, p h , t) → F u,v (h, p h , t) should be weakened, in order to cover realistic situations such as those related to cloud effects. However, it will be kept in this paper for simplicity. As shown below, this will make it possible to apply the "general" results in the previous sections. In the following sections, we will check that the theoretical analysis in Section 2 and the proposed iterative algorithms in Section 3 are valid in this framework. In particular, we consider the existence of optimal aiming strategies, their characterisation through optimality conditions and finite dimensional approximations and penalty and duality-penalty solution methods are described.

Existence and optimality
The following result holds: and Then, one has for all h and (t, p) and whence p belongs to a bounded set in H 1 (0, T ; F). This ends the proof.
From Proposition 4.1, we get the following consequences: -Theorem 2.1 can be applied and there exists at least one optimal aiming strategy for the modelled SPT plant. -Theorem 2.3 can also be applied and, assuming that at an optimalp the constraints associated to the M i are qualified, we get the necessary conditions (2.3) for some multipliers λ 1 , λ 2 ≥ 0.

Finite dimensional approximation and iterative algorithms
In practice, as in Sections 3.1 and 3.2, we must approximate the infinite dimensional problem (1.1) corresponding to (4.3)-(4.5) and replace H 1 (0, T ; F) by X N and P 0 by P 0,N . The resulting tasks are thus to solve (3.4) and/or (3.14).
Note that, in both cases, the computations of ∇J(p N ), ∇M 1 (p N ) and ∇M 2 (p N ) are needed and this requires integrals on R, with respect to (u, v) of several functions. For this reason, it is convenient to fix a set of test points (u i , v i ), 1 ≤ i ≤ I and replace these integrals by appropriate finite sums.
The total radiation reaching the receiver at time t is therefore approximated by Similarly, the deviation from the target distribution is approximated by Finally, the constraint mappings M 1 and M 2 given by (4.4)-(4.5) are approximated as follows: (4.14)

Penalisation algorithm
After these approximations, we see that, in the case of the SPT model, the objective function in (3.4) is given by: It is not difficult to compute from this identity the partial derivatives of J N µ 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 γ k used 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 γ k until 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 γ k is different for each component, that is, for each heliostat h ∈ H. Consequently, the γ k h are computed according to the rules The gradient ascent method with projection is then: 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.

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: where ψ is defined in (3.8) and M N 1 and M N 2 are 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.15 ). 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 2 and the functionality of the previous algorithms.

A numerical experiment
The behavior of the previous algorithms is illustrated for the PS10 SPT plant located in Sanlúcar la Mayor, Seville (Spain); see http://www.abengoa.com/web/en. The field of this SPT plant has 624 heliostats arranged  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 3 and 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 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.

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 8 and 9 identifies 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 Core TM i7-7700HQ CPU @ 2.80 GHz.

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            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  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.

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 2H variables 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 H increases 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.

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 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.