Access the full text.
Sign up today, get DeepDyve free for 14 days.
Loading next page...
/lp/arxiv-cornell-university/rigid-profile-input-scheduling-under-constrained-dynamics-with-a-water-3GvnyUlcGo
References
References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.
- ISSN
- 1063-6536
- eISSN
- ARCH-3343
- DOI
- 10.1109/TCST.2020.3043443
- Publisher site
- See Article on Publisher Site
Abstract
Rigid-profile input scheduling under constrained dynamics with a water network application Adair Lang, Michael Cantoni, Farhad Farokhi, Iman Shames Abstract—The motivation for this work stems from the prob- [11]. In these works, the dynamics correspond to a discrete- lem of scheduling requests for flow at supply points along an time model, or a uniformly sampled continuous-time model, automated network of open-water channels. The off-take flows with constraint satisfaction only enforced on a uniformly sam- are rigid-profile inputs to the system dynamics. In particular, pled subset of the scheduling horizon. The issue of continuous- the channel operator can only shift orders in time to satisfy time constraint satisfaction is not addressed. Specific consid- constraints on the automatic response to changes in the load. This leads to a non-convex semi-infinite programming problem, eration of rigid-profile input scheduling subject to constraints with sum-separable cost that encodes the collective sensitivity of on discrete-time dynamics is studied in [1], [12]. end users to scheduling delays. The constraints encode the linear Irrigation channels are complex physical systems with time-invariant continuous-time dynamics and limits on the state continuous-time dynamics. For channels operating under across a continuous scheduling horizon. Discretization is used to closed-loop control, linear time-invariant models are suitable arrive at a more manageable approximation of the semi-infinite program. A method for parsimoniously refining the discretization when there is ample in-channel storage [13], [14]. This moti- is applied to ensure continuous-time feasibility for solutions of vates the approach pursued below, where attention is focused the approximate problem. It is then shown how to improve cost on optimizing scheduling decisions subject to constraints, over without loss of feasibility. Supporting analysis is provided, along a continuous interval, on the evolution of a continuous-time with simulation results for a realistic irrigation channel setup to model. In particular, a non-convex semi-infinite programming illustrate the approach. formulation of the scheduling problem is considered for Index Terms—Continuous-time dynamics, load input schedul- ing, semi-infinite programming. linear time-invariant dynamics. Scheduling problems subject to constrained continuous-time dynamics are also considered in [15]–[18]. However, these works do not explicitly consider I. INTRODUCTION important aspects of computing feasible solutions for the semi- infinite program. This motivates the focus here on ensuring The problem considered in this paper is motivated by issues constraint satisfaction over the continuous scheduling horizon. that arise in the management of automated open-channel Particular emphasis is given to the construction of direct networks for rural water distribution [1]. Specifically, the discretizations with parsimonious non-uniform sampling of scheduling problem of interest pertains to the timing of rigid- time and constraints such that solutions of the corresponding profile off-take flows at supply points, given constraints on approximate problem are continuous-time feasible. This is the transient channel water-level and flow responses to load the aim of the proposed first stage of the approach. In the changes, while minimizing the collective cost of scheduling second and final stage, the cost of the first-stage feasible delays to end users. Problems of this kind may also be relevant schedule is improved by a sequential quadratic programming in other domains, e.g., energy systems and process control. approximation of the original problem. While it is possible to Some of the most widely studied scheduling problems, such extend the approach to more general dynamics, the linear time- as job and machine allocation, typically involve only static re- invariant context of this paper enables explicit characterization lationships between the variables of interest [2]–[4]. Similarly, of the state and derivatives with respect to the scheduling within the literature on irrigation networks, most scheduling variable, at any specified time for a given initial condition. studies are limited to static (or steady-state) capacity constraint In the rest the introduction, the semi-infinite program is satisfaction, with no regard for the transient behaviour [5], [6]. formulated. Aspects of the discretization-based approxima- Scheduling problems that involve constraints on the tran- tions are then introduced. To conclude, the contributions are sient behaviour of a dynamical system are considered in [7]– summarized and an outline is given for the paper. This work was supported by Rubicon Water Pty. Ltd., the Australian Re- search Council (LP130100605, LP160100666), and a McKenzie Fellowship. A. Problem formulation A. Lang is with Dept. of Electrical and Electronic Engineer- ing, The University of Melbourne, VIC 3010, Australia, and Consider a linear time-invariant system with dynamics Rubicon Water, 1 Cato St, Hawthorn East, VIC 3123, Australia. aclang@student.unimelb.edu.au M. Cantoni is with Dept. of Electrical and Electronic Engineering, The Uni- x _ (t) = Ax(t) + Bu(t) + E v (t ); x(0) = x : (1) j j j 0 versity of Melbourne, VIC 3010, Australia. cantoni@unimelb.edu.au j=1 F. Farokhi is with the Dept. of Electrical and Electronic Engineering, The University of Melbourne, VIC 3010, Australia, n n x u In this model, x(t) 2 R is the state, u(t) 2 R is the ffarokhi@unimelb.edu.au control input, and each v (t ) 2 R is a shifted load input j j I. Shames is with Dept. of Electrical and Electronic Engineering, The Uni- versity of Melbourne, VIC 3010, Australia. ishames@unimelb.edu.au at time t 2 R . Let T := [0; T ] be the scheduling horizon arXiv:2012.01626v1 [math.OC] 3 Dec 2020 2 and let N denote the integers between a > 0 and b a Problem 1 could be reformulated in terms of a switched hy- [a;b] inclusive. Given the admissible shift interval D := [ ; ] brid model of the dynamics, in which the switching instances j j for each piecewise continuous load input request correspond to the scheduling parameters in (2a). However, it appears that doing so leads to even more challenging problems. v = (t 2 [ ; T ] 7! v (t) 2 R); j 2 N ; j j j [1;m] Firstly, the set of switch states grows exponentially in the number of users, since the evolution of the channel dynamics the scheduling problem is to select the piecewise continuous depends on the combination in which requested load profiles control input u : T ! R and load input schedule shift are applied. Secondly, the natural reformulation does not variables 2 D to satisfy linear constraints on the state j j match better studied forms of switched dynamics optimiza- m u v tion problems [21]–[23]. Specifically, the switching sequence x(t; u; (v ; ) ) = x (t) + x (t ); t 2 T ; (2a) j j j j=1 0 j would not be pre-determined as in [21], this reformulation does j=1 not decompose in the way considered in [22], and it does not where seem possible to comply with the rigid profile requirements within the framework of [23], where switching sequence and x (t) := exp(At)x + exp(A(t s))Bu(s)ds; (2b) 0 switching instants are decision variables together. As such, switched system models are not considered further. Instead, x (t) := exp(A(t s))E v (s)ds: (2c) j j direct discretization approximations are pursued for obtaining a first-stage feasible solution, used to initialize a continuous The best choice, among the feasible possibilities, minimizes variable local approximation method for the original SIP. f ( ; : : : ; ) := h ( ). The differentiable h : D ! 1 m j j j j j=1 R represents the sensitivity of user j to shifting its request B. Direct discretization-based approximations (e.g., production cost impacts of delayed delivery). PROBLEM 1 (FIXED- PROFILE LOAD SCHEDULING): Given A common approach to SIPs is to relax the problem by set U of admissible control inputs, and state constraint param- sampling the constraints across the time horizon [19]. Even n n n c x c eters C 2 R and c 2 R , solve so, Problem 1 remains difficult. Hence related work has also considered discretization of the decision spaces to yield an f := min f ( ; : : : ; ) (3a) 1 m u;( ) integer linear program [1], [12]. Whilst still difficult, standard j=1 methods and solvers exist for such problems; e.g., see [24]. On s:t: Cx(t; u; (v ; ) ) c; t 2 T ; (3b) j j j=1 the other hand, shortcomings of this existing work relate to the 2 D ; j 2 N ; (3c) j j [1;m] use of uniform discretization, which can lead to unnecessarily u 2 U: (3d) large problems in terms of the number of decision variables and number of constraints, and the lack of guarantees on Several factors make it challenging to solve (3) in general: continuous-time feasibility of the solutions obtained. Aspects 1) non-convexity of the constraint set with respect to the of the discretization approach are now presented as a precursor decision variables ( ) ; j=1 to a summary of the main contributions made in this paper. 2) infinitely many constraint in (3b); and Replacing D with the finite discretized decision space 3) infinite dimensionality of the control input u. (1) (N ) m j The first factor relates to the influence of ( ) on the state D := f ; : : : ; g D (4) j=1 j j j j x(t; u; (v ; ) ) in the inequality constraints (3b); see [17]. j j j=1 for j 2 N yields the problem [1;m] Non-convexity can make it difficult to distinguish between potentially multiple local minima and the global minimum. ^ f := min f ( ; : : : ; ) (5a) 1 m ( ) With the second factor, the problem can thus be classified j=1 as a non-convex semi-infinite program (SIP); see [19] for a s:t: Cx(t; ( ) ) c; t 2 T ; (5b) j=1 comprehensive overview of such optimization problems. The 2 D for j 2 N ; (5c) j j [1;m] third factor relates to the complexity of solving constrained continuous-time optimal control problems. Given the linear where explicit dependence of x on the fixed u and given dependence of (2a) on u, this can be largely overcome requests v has been dropped for convenience. In principle, by restricting the set of admissible controls to be finite- problem (5) could be solved by exhaustive search as the dimensional, which is a commonly employed technique [20]. decision space is finite. A more sophisticated approach is Doing so yields a problem in which the main difficulty rests presented in Section II. The approach is based on iteratively solving the following related finite-dimensional optimization with how the state dynamics are affected by the decision variables ( ) ; i.e., factors 1) and 2). As such, the rest of problems. j=1 the paper is focused on these difficulties, while simplifying the The optimization problem (5) is a restriction of (3), and developments by assuming that the control input is constant; thus, f f . In fact, it is possible that (5) is infeasible for n m u ^ i.e., U = fu = (t 2 T 7! u 2 R )g. Indeed, u is instances of the finite decision sets (D ) . Under Assump- 0 j j=1 subsequently dropped as a decision variable. tion 1.1, this can be overcome by adding elements to these sets. ASSUMPTION 1.1 (S TRICT FEASIBILITY): There exists at Parsimonious augmentation is desirable as the complexity of least one schedule of load request shifts ( ) , with 2 D , solving (5) grows with the cardinality N of D , j 2 N . j j j j j [1;m] j=1 such that Cx(t; u ; (v ; ) ) c < 0 for all t 2 T . A corresponding method is proposed in Section II-B. 0 j j j=1 3 Replacing T in (5) with a finite subset of sample times at extended to accommodate combined discretization of the de- which the state constraints must hold yields a relaxation. With cision space and the constraints. The resulting discretizations are not necessarily uniform, unlike the purely discrete-time (1) (T ) T := ft ; : : : ; t g T ; i 2 N ; (6) i [1;n ] i i c formulations considered in [1], [12]. This can lead to smaller integer linear programming problems in the first-stage. In the problem becomes the second-stage, an SQP method from [26], [27] is used ;L to improve the cost without loss of feasibility. As in [17], f := min f ( ; : : : ; ) (7a) 1 m ( ) j=1 which develops a continuous-variable penalty-function based method that does not necessarily lead to improvement of s:t: C x(t; ( ) ) c ; t 2 T ; i 2 N ; i j i i [1;n ] j=1 c the cost, analytic formulations of the derivatives required (7b) to implement the SQP method are possible by virtue of 2 D ; j 2 N ; (7c) j j [1;m] the linear time-invariance of the underlying dynamics. The proposed algorithm is ultimately demonstrated on a non-trivial where C and c denote the i-th row of C and i-th entry of c, i i simulation example that is based on models of an Australian respectively. The cardinality of each constraint discretization ;L ^ ^ irrigation channel that are used operationally in the field and set T is denoted by T , which may be zero. Further, f i i historically realistic demand profiles. f , by definition. Problem (7) is tractable in the sense that it can be trans- formed into an integer linear program. This is shown in D. Outline Section II-C. However, it is possible for solutions of (7) to The rest of the paper is organized as follows. Section II in- violate (5b), and thus, be infeasible for problem (3). In view cludes discussion of the SIP discretization procedure from [25] of Assumption 1.1, and continuity of the state with respect and the modifications made to ensure finite termination when to t 2 R , this issue can be overcome by adding sample it is applied to (5) from an initially infeasible discretization of times to the sets T T and tightening the constraint (7b). the decision spaces. The continuous-variable SQP approach Specifically, the tightened problem takes the form to improving the first-stage feasible schedule is developed ;U in Section III. A formal characterization of the combined f := min f ( ; : : : ; ) (8a) 1 m ( ) j=1 two-stage algorithm is presented in Section IV. Supporting m g analytical results are provided throughout. Numerical results s:t: C x(t; ( ) ) c ; t 2 T ; i 2 N ; i j i i j=1 [1;n ] based on non-trival irrigation channel scenarios are discussed (8b) in Section V. The paper is concluded with discussion of future 2 D ; j 2 N ; (8c) j j [1;m] research directions in Section VI. g ;U ;L ^ ^ with > 0. This is a restriction of (7), whereby f f . Problem (8) is neither a restriction nor a relaxation of (3). II. S TAGE 1 - DISCRETIZATION W ITH FEASIBILITY However, as suggested above, suitable constructions of T and The Hybrid-SIP algorithm (HSIPA) from [25], which builds can be made so that solutions of (8) satisfy (5b), and thus, on [28], [29], provides a mechanism for generating constraint ;U ;L ^ ^ ^ (3b). In this case, f f f . To this end, a recent discretization sets T , and a constraint restriction parameter , algorithm for general SIPs from [25] is adapted to (5), as such that an optimizing schedule for (8) is feasible for (5) with detailed in Section II. ;U f ^ ^ f f less than a specified tolerance > 0. The HSIPA requires the following input parameters: C. Contribution initial constraint discretization sets T T , i 2 N ; i [1;n ] initial constraint restriction 2 R in (8); and >0 The main contribution is a two-stage approach to Problem 1; desired optimality tolerance 2 R for problem (5). >0 see (3). The first stage involves iterative construction of the ^ ^ discrete decision spaces D , constraint discretization sets T , Given these inputs the HSIPA iteratively determines upper and j i g ^ and constraint restriction parameter in (8), from initial lower bounds for f . These bounds converge to within the values. The resulting finite-dimensional problem yields a good specified tolerance of f . The bounds are generated by schedule of shifts for which (3b) holds; i.e., a schedule successively solving iterations of the following three subprob- lems: i) lower-bound problem (7); ii) upper-bound problem that is close in objective value to the optimal value in (5), (8); and iii) refinement problem: where only the decision space is discretized, and feasible for the original problem. The aim of the second stage is to := min (9a) improve this schedule while maintaining feasibility. To this m ( ) ; j=1 end, the original continuous decision space is re-instated, and a sequential quadratic programming (SQP) method is used to s:t: h ( ) f 0 (9b) j j approximately solve (3), initialized from the feasible schedule. j=1 The proposed two-stage approach goes beyond prior work m C x(t; ( ) ) c ; t 2 T ; i 2 N ; i j i i [1;n ] j=1 c in ensuring continuous-time feasibility via parsimoniously (9c) constructed discretizations. In particular, a method for general 2 D ; j 2 N ; and 2 R; (9d) j j [1;m] SIPs from [25] is adapted to the scheduling problem and 4 [25] This paper Description for given target objective value f > 0. The results are (SIP) (5) SIP to be solved used to update the constraint discretization sets (T ) , the f (x) f ( ; : : : ; ) Objective function i 1 m i=1 g R ^ f f Optimal objective value restriction parameter , and target objective f , as described x ( ) Decision variable(s) j=1 subsequently. Iterations proceed, until the upper and lower X (D ) Decision space f j j=1 bounds obtained lie within of each other. Each solve of g(x; y) Cx(t; ( ) ) c Inequality constraint(s) j=1 one of the three subproblems results in a candidate scheduling Y T Index set of infinite constraints m m solution ( ) . For given candidate solution ( ) , the ^ j j Y T ; i 2 N Constraint discretization set(s) j=1 j=1 i [1;n ] maximum level of constraint violation for each constraint (LBD) (7) Lower-bound sub-problem (UBD) (8) Upper-bound sub-problem i 2 N is given by: [1;n ] (RES) (9) Refinement sub-problem (LLP) (10) Lower level program m m g (( ) )) := max C x(t; ( ) ) c : (10) j i j i i j=1 j=1 t2T TABLE I A MAP OF NOTATION/LABELS IN [25] TO THE NOTATION IN THIS PAPER If g (( ) )) > 0, then the candidate schedule is not feasible i j=1 for (5). In this case, a point in time that corresponds to this maximum level of constraint violation is added to the con- straint discretization set. Specifically, given such a candidate 1) HSIPA notation and parameters: Table I relates the scheduling solution ( ) , the constraint discretization set is labels and notation used in this paper to those used in [25]. j=1 updated as follows: Notice, in particular, the multiplicity of constraint discretiza- ( tion sets here, one for each row of C in (5b), compared to m m T [ft (( ) )g if g (( ) )) > 0 i j j i j=1 i j=1 one set in the formulation of [25]. The HSIPA involves many T ; (11) T otherwise parameters. The parameters listed at start of Section II are important within the specific context of the subproblems (7) where and (8). The following parameter is introduced here to enable m m m a decision space update procedure to be activated when (5) is t (( ) )2T (( ) ) := arg max C x(t; ( ) ) c j j i j i i j=1 i j=1 j=1 t2T infeasible: (12) max k 2 N, the maximum number of consecutive con- is a corresponding maximizer of the maximum level of con- straint restriction updates allowed in the upper-bound straint violation. Feasibility of a schedule with respect to the procedure. infinite constraints (3b) corresponds to 2) Lower-bound: The lower-bound procedure of the HSIPA corresponds to Lines 2–11 of [25, Algorithm 2]. It pertains to g (( ) )) 0; i 2 N : (13) j [1;n ] i j=1 c solving the relaxed problem (7) to determine a lower-bound ;L ^ ^ The update mechanisms of the constraint restriction parameter f for f . Given a feasible solution of (7), the constraint g R and target objective f are discussed within the next discretization sets T are updated as per (11). This update, and section. The optimality tolerance should be chosen to those made in subsequent procedures of the HSIPA, ensure reflect the desired optimality gap for solving (5), which likely that successive runs of the lower-bound procedure result in a depends on the given setup. bound that converges to within a specified optimality tolerance Algorithm 2 in [25] is extended here in two ways. The of f . To ensure finite termination of the HSIPA for an initially first extension provides a way to deal with the finite sets D , infeasible problem (5), the lower-bound procedure is modified and thus, a potentially infeasible problem (5) for the given here to initially check if (7) is feasible. When it is not, the algorithm initialization. The introduction of a new parameter procedure terminates without updating the sets T , after setting enables interlacing of the HSIPA with an update procedure a flag that is used in the extension of Section II-B to trigger for the sets D that achieves eventual feasibility of (5) from j an update the discrete decision space D . an initially infeasible discretization. The corresponding update 3) Upper-bound: Lines 12–24 of [25, Algorithm 2] consti- procedure is presented in Section II-B. Secondly, the develop- tute the upper-bound procedure of the HSIPA. This procedure ment as presented here elucidates the handling of a finitely relates to solving (8), which restricts the discretized constraints indexed set of infinite constraints (i.e., n > 1). by . If (8) is infeasible, then the restriction is gradually g g g reduced, in steps =r for specified r > 1, until it becomes feasible. The modification made here, relative A. HSIPA for the scheduling problem max to [25], is to limit the number of such steps to k , Next, the notation in this paper is mapped to that used in after which the upper-bound procedure terminates and the [25], and the key algorithm parameters are identified. Aspects HSIPA returns to the lower-bound procedure, which could of the procedures associated with the aforementioned lower- itself terminate as infeasible, triggering augmentation of the bound, upper-bound and refinement problems are described decision spaces D . When (8) is feasible, and the resulting ;U in relation to generating the required constraint discretization schedule satisfies (13), the corresponding f is an upper- ^ ^ sets T and constraint restriction parameter . Generation bound for f ; i.e., the infinite constraints (5b) are feasible. of the required discrete decision spaces D is the topic of The upper-bound procedure then terminates, after a further Section II-B. Modifications made here in adapting the HSIPA step of reduction, and the HSIPA proceeds to the refinement from [25] to the scheduling problem are highlighted below. procedure. Otherwise, when (13) does not hold, the constraint 5 discretization sets T are updated, as per (11), and the upper- has an objective value less than f , the maximum constraint bound procedure is repeated, including such updates, until an restriction possible is shown as distance from c to the open (2) upper bound is found for f . These updates, and those of blue dot at t . In all cases, the problem is not SIP-feasible successive subprocedures of the HSIPA, eventually lead to an since g is greater than 0. As such, the time sample t would i i upper-bound that lies within a specified tolerance of f , as be added to the sampled time-horizon for this constraint for established in [25, Lemma 4]. the next iterate. 4) Refinement procedure: The role of the refinement pro- cedure is to improve both the upper- and lower-bounds whilst B. Extended HSIPA (EHSIPA) avoiding over-population of the sets T . The approach is Theorem 2 in [25] establishes conditions for finite termi- adapted from [30] in the HSIPA developed in [25]. The sub- nation of the HSIPA at an -optimal solution of the SIP procedure corresponds to Lines 26–42 in [25, Algorithm 2]. at hand. For the result to apply here, the corresponding A target objective value, f , is selected as the mid-point problem (5) must be feasible. This may not be the case for between the current upper and lower bounds. The correspond- certain initializations of the decision space discretizations. The ing refinement problem (9) is then solved and updates are modifications made to the HSIPA presented in the previous made on the basis of the outcome. The following lemmas section ensure finite termination, without a solution, when reveal how solving (9) may be used to improve the upper- (5) is infeasible. In this case, the HSIPA sets a flag used to ;U ;L ^ ^ or lower-bounds, f and f , respectively, as detailed in trigger augmentation of the decision space dicretizations. This the subsequent remarks. continues until feasibility of (5) is achieved . The subsequent LEMMA 2.1: If in (9) satisfies < 0, then f < f . run of the HISPA then augments the constraint discretization Proof: See Appendix A-A. sets until the sub-problems (8) yields an -optimal solution LEMMA 2.2: If an optimal schedule ( ) for (9) satisfies j j=1 of (5), as described in the preceding subsection. (13), then f f . Given decision space discretizations D = Proof: Since the resulting schedule satisfies (13), it is (1) (N ) f ; : : : ; g, N finite, j 2 N , the update j [1;m] j j feasible for (5), and hence, f f ( ; : : : ; ). By the 1 m mechanism is defined by the following: R R constraint (9b), f ( ; : : : ; ) f . As such, f f . 1 m ( ) ( ) (1) (N ) REMARK 2.1: The refinement procedure is repeated whilst + + j j j ^ ^ D D [ [ j j the conditions of Lemma 2.1 are satisfied, updating the lower- 2 2 bound to f accordingly for each such run. This can improve ( ) (l 1) (l) the overall computational time of the HSIPA as each run halves j j [ : l 2 N : (14) [2;N ] the difference between the upper and lower bound. REMARK 2.2: Satisfaction of (13) in Lemma 2.2 implies 2 R ; which is thus an upper-bound for the smallest LEMMA 2.3: Under Assumption 1.1, a finite number of constraint restriction that retains feasibility. Updating ac- decision space discretization updates (14) leads to a feasibile ;U cordingly yields improvement of f in subsequent runs of problem (5). the upper-bound procedure, as described in [25]. Proof: See Appendix A-B. REMARK 2.3: It is possible that neither Lemma 2.1 nor REM ARK 2.4: Each update (14) leads to a doubling of the Lemma 2.2 apply, leading to Line 36 of [25, Algorithm 2]. In cardinality of every D , and consequently larger integer linear this case, the constraint discretization is updated, as per (11), programs to solve in the HSIPA. It is of interest to devise a and (9) is re-solved. To moderate growth in the cardinality more parsimonious mechanism, which is the topic of future of the constraint discretization sets T , this process terminates work. Dependence of run time on the discretised decision after a specified l 2 N consecutive solves of (9), and the space cardinalities is explored empirically in Section V. max HSIPA returns to the lower-bound procedure. ASSUMPTION 2.1 (S OLUTION TOLERANCES): The feasi- In practice, the optimization problem (9) can only be solved bility of problems (7) and (8) is assessed without tolerance RES to a specified tolerance, say > 0. So as outlined in [25, to infeasibility. When found to be feasible, these problems Proposition 1] it is possible to get an inconclusive result. In are solved to within specified global optimality tolerances LBP UBP this case the HSIPA returns to the lower-bound procedure. ; 2 R , respectively. Further, it is possible to >0 5) HSIPA Illustration: Fig. 1 shows an example trajectory solve problem (10) to arbitrary optimality tolerance. of a constraint across a continuous scheduling horizon in order REM ARK 2.5: In principle, (7) and (8) can be solved to illustrate the key differences between problems (7), (8) and exactly by exhaustive search. However, this would be im- (9). The original problem (3) requires this trajectory to remain practical for any reasonably sized problem. As shown in the under the blue line (c ) for the entire time-horizon. The lower i next subsection, the discretization of decision spaces enables bound problem (7) is feasible if the constraint trajectory lies transcription to standard linear integer programs. While still below the blue dots at the specified sampled time-points, which difficult to solve, existing algorithms for such problems offer holds true in this example. For the upper-bound problem (8) much greater efficiency than exhaustive search; see [24], [31]. the constraints are still enforced only at specified sample times, Problem (10) can be solved in principle by simulation of the but the trajectory needs to be under the red-dashed line c . linear time-invariant dynamics with sufficient accuracy. (2) In this example, the constraint at t is violated. When solving With Lemma 2.3, the conditions required to apply [25, the refinement problem (9), assuming this example trajectory Theorem 2] are met. Its application yields the following result. 6 Fig. 1. An example constraint trajectory to highlight the difference between the subproblems within the HSIPA. 2 3 2 3 (1) x (t ) PROPOSITION 2.4: Given desired optimality tolerance (I C ) T 1 j1 0 i 6 7 (2) 6 7 u for problem (5), under Assumption 2.1, if there exists ~ 2 (I C ) 6 7 T 2 j2 x (t ) 6 7 0 i 6 7 f f f f m ^ = 6 7 ; b = ; j . i R such that ~ f ( ; : : : ; ) f for a schedule ( ) 6 7 >0 . 1 m j j=1 . 4 5 4 . 5 with max g (( ) )) < 0 (i.e., strictly feasible) and i2N i j=1 (T )) [1;n ] j u i c (I C ) T n jn n c c x (t ) 0 i 2 3 f f LBP UBP > ~ + + ; (15) (1 c ) (I C )b T 1 T 1 1 1 1 6 7 (1 c ) (I C )b T 2 T 2 2 2 2 6 7 then the EHSIPA terminates after finite steps at a schedule c = ; 6 7 m f 4 5 ( ) that satisfies (3b) and f ( ; : : : ; ) f + . j j=1 1 m . Proof: First note that T and D are compact sets; the (1 c ) (I C )b j T n T n n n c n c c c c h i former is a bounded real interval and the latter a finite set [32]. > (1) (2) (N ) f = ; j h ( ) h ( ) h ( ) j j j Therefore, [25, Assumption 1] is satisfied. Now from (2a), note j j j that the constraints are all continuous in the schedule variables, v(k;l) (k) (l) (l) (k) as is the objective function, since the h are differentiable. where x = x (t ) with and t as defined j j i j j i Hence for the case of a single constraint, n = 1, [25, in (4) and (6), respectively. Assumption 2] is satisfied. In the case of n > 1, the multiple Now the lower bound problem (7) can be rewritten as an constraints can be replaced with a single constraint given by: equivalent binary linear program given by max (C x(t; ( ) ) c ) 0; t 2 [0; T ]; (16) X i j i j=1 i2N > [1;n ] min f z ; (17a) (z ) j=1 j=1 which is continuous since the maximum of a set of continuous functions is continuous. Hence, under Assumption 2.1, all s:t: z c 0; (17b) j j requirements for the proof in [25] to follow are satisfied. j=1 REMARK 2.6: Proposition 2.4 only guarantees the number N > z 2 f0; 1g ; 1 z = 1; j 2 N : (17c) of steps required to reach optimality is finite. It does not j j [1;m] give any guarantees on the size of this number. The chosen Problem (8) can be reformulated similarly, except c is formed tolerance can have a significant effect on the overall number with c in place of c . Similarly, reformulation of (9) i i of steps required. Instead of letting the algorithm run until results in desired optimality guarantees are achieved, the algorithm can be stopped after at least one valid upper bound has been found. = min (18a) The schedule associated with the lowest upper-bound can then (z ) ; j=1 be returned, since this schedule is feasible, i.e., satisfies (3b). m > R The sub-optimality gap achieved with such an approach would s.t f z f 0 (18b) be highly problem and configuration dependent. j=1 z c ; (18c) j j C. Sub-problems as Integer Linear Programs j=1 For the linear dynamics considered here, the HSIPA sub- N > z 2 f0; 1g ; 1 z = 1; j 2 N ; (18d) j j [1;m] problems (7), (8), and (9), can be equivalently reformulated as binary linear programs. To this end, define the following: > > > (1 ) (1 ) : : : (1 where = . T T T 1 2 n 2 3 v(1;N ) v(1;1) v(1;2) j The reformulations above generalize the approach in [1], x x x j j j 6 v(2;1) v(2;2) v(2;N ) 7 j [12], where purely discrete-time dynamics were considered. 6 x x x 7 j j j 6 7 The generalization allows for the non-uniform sampling and = ; ji . . . 6 . 7 . . . . 4 . 5 . . . different sampling sets for each constraint generated by the v(T ;N ) v(T ;1) v(T ;2) i j i i EHSIPA. Each reformulation involves a change of variables x x x j j j 7 from to z for each user j. The corresponding relationship in the schedule ( ) . Moreover, j j j j=1 between z and is given by j j m v h i C x(t; ( ) ) = C E v (t ) C Ax (t ) i j i ` ` ` i ` j=1 ` (1) (N ) j @ : : : z = : (19) j j j j Proof: See Appendix A-C. This change of variable depends on N ; T < 1, j 2 N , j i For given feasible schedule ( ) , the quadratic program to [1;m] j j=1 i 2 N , which is not the case for the original formula- solve at each step of the SQP algorithm takes the form [1;n ] tion (3) where the corresponding sets are continuous intervals. m m X X 1 d The costs and constraints in (17) and (18) are linear in the f := min p Bp + p h ( ) + h ( ) (20a) j j j j j m p2R 2 d decision variables (z ) , except for the binary requirement j j=1 j=1 j=1 on the entries of z . Dedicated linear-integer programming m > m s:t: a (t ; ( ) ) p + b (t ; ( ) ) 0; i i j i i j j=1 j=1 solvers are readily available for such problems; e.g., see [24]. t 2 A ; i 2 N ; i i [1;n ] The difficulty of solving discrete problems in non-standard (20b) forms, such as (7), (8) and (9), is highlighted in [31]. p 2 [ ; ]; j 2 N ; (20c) j j j j The reformulations (17) and (18) can be relaxed to linear j [1;m] (q) programs by replacing the binary constraint with z 2 (0; 1), kpk : (20d) j 1 q 2 N and j 2 N . This relaxation leads to much [1;N ] [1;m] The decision variable p = p p p is the update nicer problems to solve, but at the cost of losing the rigidity 1 2 m direction from the current iterate ( ) to a new schedule, constraint, as discussed in [1]. So it is only applicable where j j=1 the load-profile is not constrained to a fixed shape. m m m a (t ; ( ) ) = r m C x(t ; ( ) ) 2 R i i j ( ) i i j j=1 j j=1 j=1 and III. STAGE 2 - L OCALIZED COST IMPROVEM ENT m m b (t ; ( ) ) = C x(t ; ( ) ) c 2 R i i j i i j i j=1 j=1 The EHSIPA described above involves discretization to ar- are constants, for the given feasible schedule and the finite car- rive at a schedule ( ) that is feasible for (3). This schedule j j=1 dinality sets A contain the maximizers T (( ) ) defined f i j i j=1 is -optimal with respect to the restriction (5) associated with in (12) for i 2 N . The positive scalar is a trust region [1;n ] ^ ^ c the final discretized decision spaces D . Since D D by j j j for the current iterate. The matrix B is a symmetric positive definition, it may be possible to improve the sum of user definite approximation of the Hessian of the Lagrangian, with sensitivities to scheduling delay, f ( ; : : : ; ), by returning 1 m jA j n i c multipliers (( ) ) , given by il l=1 i=1 to the original continuous interval decision spaces. To this end, jA j an SQP algorithm is used in a second stage of the proposed m i n L ( ) ; (( ) ) (21) j il j=1 l=1 i=1 approach to Problem 1. Starting with the feasible schedule n jA j m c i generated by the EHSIPA, the approach leads to localised X XX (l) = h ( ) + C x(t ; ( ) ) c ; j j il i j i improvement of the cost whilst maintaining feasibility with i j=1 j=1 i=1 l=1 respect to the original constraints in (3). where jA j denotes the cardinality of A , i 2 N . i i [1;n ] REM ARK 3.1: Note (20) is feasible provided the schedule A. SQP approximation around which the problem is approximated is feasible. This SQP methods involve solving quadratic programs formu- can be seen by setting p = 0. lated to approximate the original problem around each iter- The quadratic program (20) is a local approximation of the ate [33]. Adaptations of such methods to SIPs, like prob- original problem (3) relative to the schedule ( ) . The SQP j=1 lem (3), can be found in [19], [34]. The usual approach is algorithm (SQPA) proceeds by updating this schedule to the to discretize the infinite constraints. However, such relaxation schedule ( + p ) , where p solves (20) and a line search j j=1 may result in solutions that are infeasible for the original is performed to find the largest 2 [0; 1) such that problem. This could be overcome by refining the constraint g (( + p ) ) 0; 8i 2 N ; (22) j [1;n ] i j j=1 c discretization, in a similar fashion to the decision space update (14) in the EHSIPA, for example. However, this can lead to a and large number of constraints. The approach here is to sample m m m X X X constraint i 2 N at the maximizers T (( ) ) defined [1;n ] j h ( + p ) h ( ) + p h( ); (23) c i j=1 j j j j j j j m d by (12) for the given feasible schedule iterate ( ) . When j j=1 j=1 j=1 j=1 the number of maximizers is finite, as subsequently assumed, where 2 (0; 1) is a constant algorithm parameter. Corre- a tractable SQP algorithm ensues [26]. sponding updates of a , b , and A , i 2 N , h , h , i i i [1;n ] j j dt To construct the approximate quadratic program at each j 2 N , B, and are then made to become consistent [1;m] step, differentiability of the constraints is required. For the with the new iterate of the schedule. These steps repeat until problem (3) at hand, these derivatives can be formulated a stopping criterion is met. Here this criterion relates to the explicitly. step size becoming smaller than a tolerance > 0; i.e., LEMMA 3.1: The left-hand side of constraint m s C x(t; ( ) ) c 0, i 2 N , in (3) is differentiable kp k < : (24) i j i [1;n ] 1 j=1 c 8 To guarantee a unique solution to (20), the matrix B must IV. C HARACTERISTICS OF THE TWO- STAGE APPROACH remain positive definite. Since the Lagrangian (21) is non- The two-stage approach to Problem 1 is as follows. First, convex in the schedule a standard Broyden-Fletcher-Goldfarb- the EHSIPA described in Section II is applied to generate Shanno (BFGS) update formula may lead to non-positive f a schedule that is feasible for (3) and -optimal for the definite matrices. Here a damped BFGS update formula from restriction associated with a discretization of the decision [33, p.537] is used. The damping ensures a curvature condition spaces. In the second stage, this SIP feasible schedule is then is satisfied and that B remains positive definite if the initializa- improved in cost by application of the SQPA described in tion is positive definite. The update formula from [33, p.537] Section III. requires an estimate of the Lagrangian multipliers at current THEOREM 4.1: Under Assumptions 1.1 and 2.1, suppose f LBD UBP iterate which can be obtained in the process of solving (20). ; ; and satisfy (15). Then the two-stage approach The trust region is used to reflect confidence in the terminates finitely with a schedule ( ~ ) such that i) ~ 2 j j=1 j previous approximation. This can be dynamically updated to D , j 2 N , ii) (3b) holds, and iii) j [1;m] allow for larger steps when the approximation is deemed more f f f ( ~ ; : : : ; ~ ) f + ; (26) accurate, and refined as the approximation becomes coarse, in 1 m order to reduce the number of iterations needed to perform the where f and f are defined by (3) and (5), respectively, and line search. Here, [33, Algorithm 4.1] is used to update the f := min f ( ; : : : ; ); s.t 2 D ; j 2 N : (27) 1 m j j trust region. [1;m] ( ) j=1 THEOREM 3.2: Given > 0 and 2 (0; 1), suppose the Proof: Proposition 2.4 gives the right hand side of (26). SQPA is initialized with > 0, B positive definite, and a Theorem 3.2 ensures that the SQPA maintains or improves the schedule that is feasible for (3) with cost f . Then cost whilst maintaining SIP feasibility. Hence, the bound and i) every schedule iterate is feasible for (3), and feasibility still hold after stage 2. The left-hand side of (26) ii) the algorithm terminates in a finite number of steps, such follows because the cost of any feasible schedule is an upper that when the initial solution of (20) is non-zero, bound for the optimal cost and that the optimization problem defining f is a relaxation of (3); i.e., there are no constraints h ( ) < f (25) j j on the dynamics. j=1 REM ARK 4.1: Theorem 4.1 holds for any second stage at termination. algorithm that takes an initial schedule that is feasible for (3) and terminates after a finite number of steps with a schedule Proof: See Appendix A-D satisfying (25) and (13). REMARK 3.2: Either the SQPA terminates immediately, REM ARK 4.2: The inequality (26) bounds the ob- with the initial schedule, or a new SIP feasible schedule with tained optimality gap. Specifically, f ( ; : : : ; ) f 1 m strictly improved cost is found in finite steps. By contrast, strict f ( ; : : : ; ) f f + f . If an upper bound is 1 m improvement cannot be guaranteed by the penalty methods known for f then the last expression can be calculated a proposed in [17]. priori. However, in practice these bounds are conservative REMARK 3.3: If the requested load inputs v are restricted since f f . This bound relies on calculation of (27) which to continuous functions, and certain regularity assumptions should be easier than (3) since there are no constraints from hold, and the sets A at each step are modified to include the dynamics. In particular if h are convex then (27) is a points that can be extended (in an appropriate sense) to each j convex problem which could be solved using standard convex of the global maximizers of another schedule that is in the optimization methods such as those in [35]. neighborhood of the iterate ( ) , then it is possible to show j=1 REM ARK 4.3: The main result relies on Assumption 2.1. If that the SQPA, with sufficiently small > 0, converges to there exist algorithms that satisfied Assumption 2.1 with the a point that satisfies the first order optimality conditions of original continuous decision spaces, then the HSIPA algorithm (3); i.e., to a local-stationary point. This result is established for constraint discretization would find the optimal solution, in [26]. However, improvement in the cost is the primary aim negating the need for the second localized SQPA based re- here, and since it is difficult to show the regularity assumption finement stage. The decision spaces are discretized to arrive at holds and to ensure A contains the necessary points, the tractable problems that satisfy Assumption 2.1. The localized conditions are relaxed to those in Theorem 3.2, at the expense refinement stage can then lead to improvement. of not ensuring local optimality. Although Theorem 4.1 guarantees finite termination it does REMARK 3.4: The algorithm presented in [26] also allows not say anything about the convergence rate or computa- for infeasible initial schedules via an exact penalty method. tional complexity and specifically how to choose the initial However, it is observed in practice that this approach is very discretizations of the decision spaces and constraints. The sensitive to initialization. Within the context of the numerical subsequent discussion relates to why a good choice of initial examples presented in Section V, initialization of the approach discretizations is not obvious. at realistic original flow requests failed to yield a schedule that satisfies (13). This is not unexpected since stationary A. Effect of Initial Discretizations points for the exact penalty reformulation of the problem are not necessarily feasible. This motivated development of the The two-stage approach based on EHSIPA and SQPA is EHSIPA. guaranteed to converge for any choice of initial discretiza- 9 tions of the decision spaces and constraints in (3). However, mentioned complexities regarding the initial discretizations this choice can have significant impact on both the overall are demonstrated for a practical example from automated computational time and the achieved optimality gap. irrigation networks in the subsquent section. 1) Discretization and Computation Time: The integer pro- grams in the EHSIPA have dimensions linked explicitly to N V. NUMERICAL RESULTS and T , the cardinalities of the discretized decision spaces and constraint discretization sets. Hence, it may seem intuitive to In this section, a realistic scheduling setup is investigated start with these being small in order to reduce the solve time numerically based on operational data for an irrigation channel of these sub-problems. However, sparse initial sets may lead in south-eastern Australia. Results obtained from applying the to slower overall computation time. first stage of the approach developed above are presented To see this, first consider sparse initial constraint discretiza- for nine instances of the initial discreizations of the decision tion sets T . Using these may lead to an initial lower-bound space and constraints. These initializations range from fine within the EHSIPA that is far from the optimal, e.g., if the uniform discretizations of both the decision space and the constraint set were empty, then the initial lower bound would constraints to coarse discretizations of both. The outcomes are be f . Further, the first call to the upper-bound procedure may compared in terms of feasibility, optimality and computational require many iterations to add sufficient points to the constraint burden. The results reveal that a dense initialization does not necessarily yield the best outcome. For the example considered discretization sets to find a feasible schedule. it is observed that coarse initialization can yield a feasible Second, consider sparse initial discretized decision spaces schedule that is as good as one obtained from a much denser D . This could yield an infeasible problem (5) and several starting point, without the optimization sub-problems involved iterations of the EHSIPA could be required before (5) becomes becoming overly large. The effectiveness of the second stage feasible. Additionally, a sparse decision space may increase the is also explored, including comparisons with the penalty based likelihood that the resulting schedule from the first stage is far gradient methods from [17], in terms of improvement in cost from a stationary point of the original problem (3). This can and computational burden. increase the convergence time of the SQPA. These issues are discussed further within the context of a particular example in the numerical results Section V. A. Application Gravity Fed Irrigation Networks 2) Discretization and Optimality Gap: The SQPA is guar- anteed to improve the initial schedule if possible. However, An irrigation channel is a cascade of pools. Following [13], due to non-convexity of (3), the iterates tend towards local [36], each pool is modeled in the Laplace domain as stationary points of (3). In general, there is no measure of c c c in;i out;i out;i t s d;i y (s) = e q (s) q (s) o (s); i i i+1 i how far these local stationary points are from the global s s s optimum. Hence, the sub-optimality gap achieved for the two- p where c and c (in 1=(min m)) are discharge rates in;i out;i stage approach is strongly linked to the sub-optimality of the determined by the physical characteristics of the gates used to first stage solution; i.e., the distance between f and f . set the flow between neighbouring pools, and t (in min) is d;i REMARK 4.4: Nothing can be said about the achieved the delay associated with the transport of water along the pool. optimality gap if the SQPA is initialized at an arbitrary feasible Here, o (s) denotes the overall off-take load on pool i, that is, schedule, rather than an optimizer for (5). the sum of all the water supplied to the farms connected to this The difference f f depends on the discretization of the pool. Moreover, q (s) represents the head over the upstream decision spaces. This difference can be made arbitrarily small gate of pool i raised to the power of 3=2 which is proportional by adding sufficiently enough points to the discretzized de- to the flow (in m =min) of water from pool i 1 to pool i and cision spaces (D ) . However, this approach is impractical j=1 y (s) denotes the water level (in m) in pool i. For the purpose since this may lead to prohibitively large integer programs to of this example, the delays are replaced with a first-order Pade ´ be solved as part of the HSIPA. This is highlighted further in approximation . Each pool is controlled, locally, by Section V. As mentioned in Remark 2.4, it is of interest to investigate ( s + 1) i i q (s) = (u (s) y (s)) + q (s); i i i i i+1 the use of additional information to inform the discretization(s) s( s + 1) of smaller decision spaces. Note that it is possible for a where , , and are appropriately selected control smaller discretization to result in a smaller optimality gap. For i i i i parameters. Furthermore, u (s) (in m) denotes the water- instance, consider the case where the discretization for each i level reference signal of pool i. The model for a series of user is chosen with cardinality equal to one, containing only ^ pools using the aforementioned model, with first-order Pade ´ the shift from an optimal schedule for (3); i.e., D = f g approximation, can be represented in the form of (1) where for j 2 N . In this case, the smallest optimality gap is [1;m] there are four states per pool, two for the controller and two achieved even though the discretizations are very sparse. for the level dynamics. A possibility in this direction is to re-run the EHSIPA after the first feasible schedule is found, with a finer discretization Note that the choice of a first-order Pade ´ approximation is justifiable as the centered at the solution of the previous iterate such that pool delays are all parts of closed-loops (with local controllers), with loop- the cardinality of the decision spaces remain the same. This gain cross-overs that are sufficiently small to make the overall closed-loop preliminary idea, along with the effects of the other afore- behavior insensitive to the approximation error [14]. 10 TABLE II TABLE III PARAMETERS, c , c (IN 1/(min m)), t ( IN min), AND UNITLESS PARAM ETERS FOR OFF-TAKE LOADS USED IN THE SIM ULATION. s ; d in;i out;i d;i ij ij CONTROL PARAMETERS ; AND FOR i 2 N FOR THE (IN min) AND m MULTIPLIED BY A DISCHARGE COEFFICIENT HAS UNITS i i i [1;10] ij DYNAMICAL SYSTEM USED IN THE SIMULATION. OF m =min Pool No c c t in;i out;i d;i i i i Pool No s d m s d m i1 i1 i1 i2 i2 i2 i = 1 0:1079 0:1080 1 0:0156 46:648 3:452 i = 1 200 360 0:0645 500 360 0:0322 i = 2 0:0777 0:0777 1:67 0:0091 72:403 5:213 i = 2 200 180 0:0322 500 180 0:0322 i = 3 0:0586 0:0586 2:33 0:0065 99:274 7:084 i = 3 200 360 0:0322 500 360 0:0645 i = 4 0:1269 0:1269 1:67 0:0084 60:305 3:972 i = 4 200 180 0:0645 500 180 0:0322 i = 5 0:0313 0:0313 1:83 0:0092 110:85 8:878 i = 5 200 360 0:0322 500 360 0:0322 i = 6 0:0456 0:0507 4 0:0036 152:73 10:36 i = 6 200 180 0:0290 500 180 0:0580 i = 7 0:0725 0:0725 1:33 0:0119 64:978 4:885 i = 7 200 360 0:0580 500 360 0:0290 i = 8 0:0216 0:0216 3:67 0:0080 147:65 10:28 i = 8 200 180 0:0322 500 180 0:0322 i = 9 0:0366 0:0366 1:67 0:0100 98:231 7:816 i = 9 200 360 0:0322 500 360 0:0645 i = 10 0:2062 0:2331 2 0:0100 48:156 2:101 i = 10 800 360 0:0285 500 180 0:0285 2) Algorithm parameters for the EHSIPA and SQPA stages: In addition to the parameters mentioned throughout Section II, B. Example Parameters and Setup the HSIPA requires specification of the tolerance for problem 1) Problem Data: Table II shows the parameters for the 10 RES LLP (9), 2 R , an initial tolerance for (10), 2 R , >0 >0 pool system used in this example, which are from validated LLP which is refined throughout the HSIPA by steps models of a real channel provided by Rubicon Water Pty Ltd. LLP LLP LLP =r for specified r > 1, as detailed in [25, Algo- A fixed = 0:7 is used for all pools i 2 N . The reference i [1;10] rithm 2]. Each parameter is independently varied to determine input u is a set-point to a lower-level feedback controller. In a set of parameters that give reasonable and reliable computa- this context it is common that the reference be kept constant tional performance. This results in the set of parameters chosen as it is desired to maintain a constant level of supply to the f g g LLP LBP UBP RES LLP max as ( ; ; r ; r ; ; ; ; ; l ; k ) = max offtakes. For this setup it is fixed to 1m for all pools, i.e., 3 4 (5; 10 ; 1:5; 1:2; 0:5; 0:5; 10 ; 0:01; 10; 10). The focus of u = 1 . 0 10 the simulation study is directed to the effect of the initial The state constraints are envelope constraints on the water discretization of both the decision spaces and constraints. level such that y (t) 2 [y ; y ] for all t 2 [0; T ] where y = 0:9 i i The parameters for the SQPA are chosen similarly, resulting (m) and y = 1:075(m) 8i 2 N [f9; 10g and y = 0:88 [1;4] i in parameters ( ; ; B; ) = (5; 0:33; I ; 0:5). step m (m) and y = 1:10 8i 2 N . [5;8] i ^ 3) EHSIPA Initializations: The discretizations D are cho- ij In this example, there are two users for each pool. Each sen to be a uniformly sampled version of the continuous user j, in each pool i, has a requested off-take profile defined sets D = [ 180; 180] with sample period , i.e., D = ij ij by a start time, s , duration d (both in min), and magnitude ij ij f ; + ; : : : ; g for all i 2 N ; j 2 f1; 2g, and the ij ij ij [1;10] m , which is proportional, via a discharge coefficient, to flow ij initial discrete constraint sets T are chosen to be uniformly 3 i (in m =min). This gives typical pulse shape for off-take flows sampled time points with a period . Nine different simulation in irrigation networks, i.e., v (t) = m (H (t s ) H (t ij ij ij study configurations are considered as outlined below: s d )) where H : R ! f0; 1g denotes the continuous time ij ij 1) ( ; ) = (60; 15); Heaviside step function. Table III shows the particular off- 2) ( ; ) = (30; 15); take load parameters used in this example. The top of Fig. 2. 3) ( ; ) = (15; 15); displays the requested orders, which follow a realistic demand 4) ( ; ) = (5; 15); pattern based on historic orders. It is of note the initial requests 5) ( ; ) = (15;1) (i.e., initial constraint set is empty); cause significant violation of the constraints; see top of Fig. 3. 6) ( ; ) = (15; 30); The bounds on the allowable shifts are set to = 180 ij 7) ( ; ) = (15; 1); and = 180 for all i 2 N ; j 2 f1; 2g, i.e., the orders ij [1;10] 8) First the EHSIPA is run with ( ; ) = (30; 15) then the can be scheduled by shifting the requested load forward or EHSIPA is rerun with discretization centered around the backwards by up to three hours. The end-user sensitivity is 2 optimal solution with ( ; ) = (2:5; 15). Specifically, measured with a quadratic function h ( ) = 0:01 for all for user j on pool i denote as the optimal shift from ij users, i.e., each user experiences a greater cost the further configuration 2. The restricted set is chosen as D = away from the requested start-time. The ideal schedule for ij f 30; 27:5; : : : ; ; : : : ; + 30g; each user is to not have their order shifted at all. The planning ij ij ij ij 9) ( ; ) = (15; 0:25) with a single iteration of (7) only. horizon is set to T = 1440min (1 day), for all simulations. In open-channel irrigation networks, it is desirable to reduce The first four configurations are used to explore the effect the required lead-time for users to request desired off-take of the decision space discretizations, while configurations profiles. In practice the lead time is typically in the order of 3 and 5-7 explore the effect of the initial constraint dis- days, but there are network operators striving for lead times of cretization. Configuration 8 examines the potential for an around two hours. Hence, a computational time in the order improved method of decision space updates as discussed in of an hour for solving the scheduling problem is tolerable. Section IV-A2. Configuration 9, which does not involve update The aforementioned setup gives all the problem data needed of the initial discretizations, is equivalent to the discrete-time to formulate (3). system methods in [1]. It is included for comparison with the 11 approach proposed here. C. Implementation of proposed two-stage approach 1) EHSIPA implementation: Numerical solver ode45 in Matlab is used to evaluate the integrals in the components x (t), x (t), i 2 N ; j 2 N required to compute 0 [1;10] [1;2] ij (10). By exploiting linearity and time-invariance of (1), each constraint k is calculated for a given schedule by summing x (t) with shifted versions of x (t) and multiplying the result by the corresponding row C . The integrals are evaluated on a uniform discretized set with an initial sample period chosen as = 0:1. If, during the HSIPA, the maximum difference between any two consecutive samples of the constraint is LLP greater than the current required tolerance , the sample period is reduced accordingly and the intergrals are re- evaluated on discretization set with the smaller sample period. The sub-problems in the HSIPA are formed as the linear (mixed-)integer programs described in Section II-C and solved using Gurobi [37], interfaced with Matlab. 2) Localized refinement stage: For the second localized refinement stage the following three algorithms are compared: Fig. 2. Top: Unscheduled 20 off-take profiles. Bottom: Off-take profiles under i) the proposed SQPA, described in Section III; a feasible schedule using proposed approach with configuration 4 followed by SQPA. This feasible sub-optimal schedule is close to original, as per users ii) the projected gradient algorithm in [17] with log-penalty desire. and = 1500; iii) the projected gradient algorithm in [17] with exponential- algorithm initialized according to configuration 4. The corre- penalty and # = 1500. sponding water level trajectories are displayed in Fig. 3. The The line search in the projected gradient algorithm in [17] shifted off-takes remain “close” to the requested ones under used for ii) and iii) is modified slightly to include a feasibility the sub-optimal schedule, and the levels remain within the condition similar to (22). This allows the comparison to focus constraints. Achieving constraint satisfaction relies on the first on the final objective value and computation time. The deriva- stage, EHSIPA, to find a feasible schedule. From Fig. 4 it tives required to implement the SQPA can be constructed via can be seen that the purely discrete-time (uniform sampling) Lemma 3.1 from simulations of the linear system dynamics, approach from [1] is not able to meet the constraints even as used to evaluate (10) in the EHSIPA. By contrast, the with fine discretization (i.e., = 0:25, configuration 9). By derivatives of the augmented penalty terms needed for the contrast, the EHSIPA is able to yield a feasible solution with penalty based algorithms are computed using a first order only 95 constraints (configuration 5); see Fig. 5. In addition, in backward finite difference method; i.e., the partial derivative configurations 3, 5 and 6, which have the same initial decision of given function l() with respect to shift at the current m spaces, the proposed algorithm is almost a whole order of schedule ( ) is approximated as j=1 magnitude faster than in configuration 9; see Fig. 6. l(( ) ) l(( ) ; ) @l j j j6=` ` s j=1 The runtime of the EHSIPA increases in proportion to the ; (28) size of the decision spaces as highlighted by configurations 1- ` s 4 in Fig. 6. Configurations 3 and 5-7, have N = 24 possible with set to be smaller than the final sample period from choices for each of the 20 off-takes, which corresponds to all simulations required for the EHSIPA (for this example 24 possible combinations. As such, it is intractable to find = 0:01). This approximation is computationally much more the optimal solution via a brute-force (exhaustive search) efficient and numerically robust, compared to calculating the approach. This highlights the importance of formulating the derivatives of the penalty functions on the basis of explicit sub-problems as (mixed-) linear-integer programs as such formulae, as these do not directly correspond to a simulation problems are amendable to widely available dedicated solvers. of the linear dynamics. Configurations with the same initial decision spaces, i.e., 3 and All simulations were carried out with Matlab 2018b on a 5-7, are used to explore the effect of the initial constraint dis- Windows PC with 16GB RAM and Intel Core i7-4790K CPU cretization. Of these, configuration 6, with = 30, resulted in @4.00GHz processor. the fastest run time of the first stage (EHSIPA). Configuration 5 corresponds to the coarsest initial and final discretization, but D. Results more constraints are added during the algorithm execution. As The top of Fig 2 shows the requested (unscheduled) off-take such, it requires many more calls to the lower-level problem profiles and the bottom shows the off-take profiles shifted with (10) and iterations of the sub-problems. Note in particular, that a feasible sub-optimal schedule obtained from the proposed the initial lower-bound is much worse for configurations 5 and 12 Fig. 3. Levels y (t) (in m) in response to nominal off-takes, where constraints are clearly violated on both the upper and lower bounds (top figure); and feasible schedule from configuration 4, where the levels are within the constraint limits (bottom figure). 100 5 Fig. 4. Levels y (t) in response to the off-take schedule using the uniform discretization method from [1], i.e., configuration 9, (left), which does not meet constraints, and a off-tale schedule from configuration 3 combined with 2 SQPA, where the levels are within the constraint limits. 1 2 3 4 5 6 7 8 9 6, as shown in Fig. 7. This accentuates the points discussed Configuration No. in Section IV-A1. The configurations with finer discretization resulted in Fig. 5. Number of constraints added throughout the EHSIPA algorithm and smaller objective f obtained from EHSIPA; see Fig. 9. the total number of constraints for the 9 different configurations. Configuration 8 results in a final objective that is 0:67% lower than configuration 4, however it achieves this in 8:4% of the high fidelity models which include, low-level control system time. This highlights the potential for further improvements to actuator saturation and the St Venant PDE model for the open- decision space update methods. water dynamics [38]. Fig. 8. shows the optimal schedule only The arrows in Fig. 9 show the impact of the second stage on slightly violates the constraints for only two of the pools. the overall objective and computation time. Both a longer run time and larger improvement in objective is made with stage VI. CONCLUSIONS AND FUTURE WORK 2 for configuration 1. The logarithmic penalty method was the fastest for all configurations but only resulted in marginal A two-stage approach to finding a sub-optimal feasible improvement of the objective in each case. The proposed solution to a rigid-profile input scheduling problem for a method resulted in improvement in all configurations and continuous-time linear time-invariant dynamical system is had the greatest improvement for configurations 3, 4 and 8. proposed. The first stage pertains to the discretization of both Interestingly, for configuration 4 and 8, the exponential penalty the decision spaces and constraints. Through iteratively refin- method resulted in an increase in the overall objective, shown ing the discretizations and solving tractable (mixed-)integer by the positive gradient of the dotted arrow. That is, the penalty linear programs as sub-problems, a schedule that is feasible method cannot guarantee a strict improvement of the schedule, for the original continuous-time problem is generated. This unlike the SQPA method; see (25). schedule is then improved by the second stage in a way Finally, to highlight the applicability of the proposed method that guarantees improvement in the objective and feasibility in practice, the water levels for the unscheduled and a feasi- with respect to the original problem. A sequential quadratic ble sub-optimal scheduled off-take loads are simulated using programming method is applied for this stage and shown to # of constraints added in EHSIPA # of constraints at completion of EHSIPA 13 [4] A. Allahverdi, “The third comprehensive survey on scheduling problems (7) with setup times/costs,” European Journal of Operational Research, (8) 8 vol. 246, no. 2, pp. 345–378, 2015. [5] S. Hong, P.-O. Malaterre, G. Belaud, and C. Dejean, “Optimization (9) of irrigation scheduling for complex water distribution using mixed Other 6 integer quadratic programming (fMIQPg),” in Proceedings of the 10th International Conference on Hydroinformatics (HIC 2012), 2012. [6] J. M. Reddy, B. Wilamowski, and F. Cassel-Sharmasarkar, “Optimal 4 scheduling of irrigation for lateral canals,” ICID Journal on Irrigation and Drainage, vol. 48, no. 3, pp. 1–12, 1999. [7] M. Zachar and P. Daoutidis, “Scheduling and supervisory control for cost effective load shaping of microgrids with flexible demands,” Journal of Process Control, vol. 74, pp. 202 – 214, 2019. Efficient energy management. 1 2 3 4 [8] L. S. Dias, R. C. Pattison, C. Tsay, M. Baldea, and M. G. Ierapetritou, 10 10 10 10 “A simulation-based optimization framework for integrating scheduling Computation Time (s) and model predictive control, and its application to air separation units,” Computers & Chemical Engineering, vol. 113, pp. 139 – 151, 2018. [9] A. Allman, M. J. Palys, and P. Daoutidis, “Scheduling-informed opti- Fig. 6. The breakdown of the total runtime for the EHSIPA for each mal design of systems with time-varying operation: A wind-powered configuration. ammonia case study,” AIChE Journal, vol. 65, no. 7, p. e16434, 2019. [10] J. I. Otashu and M. Baldea, “Scheduling chemical processes for fre- quency regulation,” Applied Energy, vol. 260, p. 114125, 2020. [11] B. Burnak, N. A. Diangelakis, J. Katz, and E. N. Pistikopoulos, “In- tegrated process design, scheduling, and control using multiparametric programming,” Computers & Chemical Engineering, vol. 125, pp. 164 – 184, 2019. [12] F. Farokhi, M. Cantoni, and I. Shames, “A game-theoretic approach to distributed scheduling of rigid demands on dynamical dystems,” in Proceedings of the Australian Control Conference, pp. 147–152, 2016. [13] E. Weyer, “System identification of an open water channel,” Control Engineering Practice, vol. 9, no. 12, pp. 1289 – 1299, 2001. [14] M. Cantoni, E. Weyer, Y. Li, S. K. Ooi, I. Mareels, and M. Ryan, “Control of Large-Scale Irrigation Networks,” Proceedings of the IEEE, vol. 95, no. 1, pp. 75–91, 2007. [15] A. Flores-Tlacuahuac and I. E. Grossmann, “Simultaneous cyclic scheduling and control of a multiproduct cstr,” Industrial & Engineering Chemistry Research, vol. 45, no. 20, pp. 6698–6712, 2006. [16] J. Zhuge and M. G. Ierapetritou, “Integration of scheduling and control with closed loop implementation,” Industrial & Engineering Chemistry Research, vol. 51, no. 25, pp. 8550–8565, 2012. [17] F. Farokhi, M. Cantoni, and I. Shames, “Scheduling rigid demands on Fig. 7. Figure showing convergence of upper and lower bounds of continuous-time linear shift-invariant systems,” in Proceedings of the EHSIPA, for different initial discretizations of T . A less dense discretization, IEEE Conference on Decision and Control, pp. 5358–5363, 2015. = 1, results in more a conservative initial lower bound. Whereas, a dense [18] C. Tsay, A. Kumar, J. Flores-Cerrillo, and M. Baldea, “Optimal demand discretization, = 1 takes much longer to calculate initial lower bound. response scheduling of an industrial air separation unit using data-driven dynamic models,” Computers & Chemical Engineering, vol. 126, pp. 22 meet both of these requirements. Simulation results from a – 34, 2019. realistic example from automated irrigation networks illustrate [19] R. Hettich and K. O. Kortanek, “Semi-infinite programming: theory, the advantages and some of the properties, trade-offs for methods, and applications,” SIAM review, vol. 35, no. 3, pp. 380–429, the proposed algorithm(s). Future work is to explore how [20] J. T. Betts, Practical Methods for Optimal Control and Estimation prior information or multiple iterations of the first stage of Using Nonlinear Programming. Society for Industrial and Applied the algorithm can be used to inform more suitable choices Mathematics, second ed., 2010. [21] X. Xu and P. J. Antsaklis, “Optimal control of switched systems based of restricted decision spaces. Other work can be done on on parameterization of the switching instants,” IEEE Transactions on characterizing the convergence for the various algorithms, and Automatic Control, vol. 49, pp. 2–16, jan 2004. development of feasible methods for computation of useful [22] F. Zhu and P. J. Antsaklis, “Optimal control of hybrid switched systems: A brief survey,” Discrete Event Dynamic Systems, vol. 25, no. 3, optimality gaps to detect acceptable solutions. Finally, further pp. 345–364, 2015. work could explore how to update the schedule, perhaps on [23] T. Das and R. Mukherjee, “Optimally switched linear systems,” Auto- a receding horizon basis, to accommodate for model mis- matica, vol. 44, no. 5, pp. 1437–1441, 2008. [24] L. A. Wolsey, Integer Programming. Wiley Series in Discrete Mathe- matches and to also account for changes to requested orders. matics and Optimization, Wiley, 1998. [25] H. Djelassi and A. Mitsos, “A hybrid discretization algorithm with guaranteed feasibility for the global solution of semi-infinite programs,” REFERENCES Journal of Global Optimization, vol. 68, no. 2, pp. 227–253, 2017. [26] C. J. Price and I. D. Coope, “An exact penalty function algorithm [1] J. Alende, Y. Li, and M. Cantoni, “A f0, 1g linear program for fixed- for semi-infinite programmes.,” BIT. Numerical Mathematics, vol. 30, profile load scheduling and demand management in automated irrigation p. 723, jan 1990. channels,” in Proceedings of the 48th IEEE Conference on Decision and [27] C. J. Price, “Non-linear semi-infinite programming,” 1992. Control, pp. 597–602, 2009. [28] A. Mitsos, “Global optimization of semi-infinite programs via restriction [2] M. L. Pinedo, Scheduling: theory, algorithms, and systems. Springer, of the right-hand side,” Optimization, vol. 60, no. 10-11, pp. 1291–1308, 2016. 2011. [3] J. Blazewicz, M. Drabowski, and J. Weglarz, “Scheduling multiprocessor [29] A. Mitsos and A. Tsoukalas, “Global optimization of generalized semi- tasks to minimize schedule length,” IEEE Transactions on Computers, infinite programs via restriction of the right hand side,” Journal of Global vol. 100, no. 5, pp. 389–393, 1986. Optimization, vol. 61, no. 1, 2015. Configuration No. 14 Fig. 8. High fidelity St Venant PDE models based simulation of levels y (t) in response to nominal off-takes, where constraints are clearly violated on both the upper and lower bounds (top figure); and a feasible off-take schedule from configuration 4, where the levels are almost all within the constraint limits (bottom figure). Objective Value Vs Computation Time 1) - EHSIPA 1) - Stage 2 Log 1) - Stage 2 Exp 1) - Stage 2 SQP 2) - EHSIPA 2) - Stage 2 Log 2) - Stage 2 Exp 2) - Stage 2 SQP 3) - EHSIPA 3) - Stage 2 Log 3) - Stage 2 Exp 600 3) - Stage 2 SQP 4) - EHSIPA 4) - Stage 2 Log 4) - Stage 2 Exp 4) - Stage 2 SQP 8) - EHSIPA 8) - Stage 2 Log 8) - Stage 2 Exp 8) - Stage 2 SQP Minimum Objective Found 0 100 400 600 800 1;000 1;200 1;400 1:26 1:28 1:3 1:32 Total computational time (s) 10 Fig. 9. Runtime vs objective for the different choices of D and different algorithms for stage 2 of proposed approach ij [30] A. Tsoukalas and B. Rustem, “A feasible point adaptation of the engineering for irrigation systems: Successes and challenges,” Annual Blankenship and Falk algorithm for semi-infinite programming,” Op- Reviews in Control, vol. 29, no. 2, pp. 191–204, 2005. timization Letters, vol. 5, no. 4, pp. 705–716, 2011. [37] I. Gurobi Optimization, “Gurobi Optimizer Reference Manual,” 2016. [38] M. H. Chaudhry, Open-channel flow. Springer Science & Business [31] P. Belotti, C. Kirches, S. Leyffer, J. Linderoth, J. Luedtke, and A. Ma- Media, 2007. hajan, “Mixed-integer nonlinear optimization,” Acta Numerica, vol. 22, pp. 1–131, 2013. [32] W. Rudin, Functional Analysis. International series in pure and applied APPENDIX A mathematics, McGraw-Hill, 1991. PROOFS OF L EM MAS AND T HEOREMS [33] J. Nocedal and S. J. Wright, Numerical Optimization. Springer Series in Operations Research and Financial Engineering, New York, NY : A. Proof of Lemma 2.1 Springer New York, 2006., 2006. Assume < 0. Since, the objective of (9) is then [34] R. Reemtsen and J.-J. Ruckmann, ¨ Semi-Infinite Programming. [elec- tronic resource]. Nonconvex Optimization and Its Applications: 25, optimality of < 0 implies that for = 0 there are no m m Springer US, 1998. feasible schedules ( ) for (9). Let R := f( ) : 2 j j j j=1 j=1 [35] S. P. Boyd and L. Vandenberghe, Convex Optimization. Cambridge D ; 8j 2 N ; h ( ) f g, which is non-empty j [1;m] j j University Press, 2004. j=1 R ;L [36] I. Mareels, E. Weyer, S. K. Ooi, M. Cantoni, Y. Li, and G. Nair, “Systems since f is chosen to be greater than a lower-bound f as Objective Value 15 per [25, Algorithm 2]. Then, since = 0 is not feasible for Note that [k + 1] 2 D since (20c) must hold and through j j (9), for all ( ) 2 R there exists a constraint i 2 N the line search condition (22) then j [1;n ] j=1 c and t 2 T such that i m m g (( [k] + p [k]) ) = g (( [k + 1]) ) 0: j j j i j=1 i j=1 C x(t; ( ) ) c > 0: (29) i j i j=1 m m Hence ( [k + 1]) satisfies (13). As ( [0]) is assumed j j j=1 j=1 to be feasible, statement one follows by induction. Therefore, since T T , there is no schedule in R that i i ii) Since B[0] is postive definite, B[k] is positive definite satisfies (5b). Hence f < f which completes the proof. for all k 2 N ; see [33, Chapter 18]. Let m m B. Proof of Lemma 2.3 X X 1 d m (p) := p B[k]p + p h ( [k]) + h ( [k]): c m k j j j j j Let ( ) be a schedule satisfying Assumption 1.1; i.e., j j=1 2 d j=1 j=1 a strictly feasible schedule. Then there exists > 0 such that Since B[k] is positive definite, for p 6= 0, c m g (( ) ) : (30) i j j=1 m m X X In view (2b) and (2c), the dependence of the left-hand side m (p) > p h ( [k]) + h ( [k]): (32) k j j j j j m d of the constraint (3b) is continuous in ( ) , whereby g is j j=1 j=1 j=1 i continuous as the max of continuous functions. Therefore, for The optimal objective value for (20) is given by every > 0 there exists () > 0 such that f [k] = m (p [k]): (33) c k maxk k < () (31a) An upper bound for f [k] is m (0), because p = 0 is a feasible m c m ) g (( ) ) < g (( ) ) + + : (31b) i j=1 i j j=1 solution for (20). Combining this with (33) and (32) gives m m In particular, with = , if given schedule ( ) satisfies j=1 X X p [k] h ( [k]) + h ( [k]) < f [k] the corresponding condition (31a), then it is feasible for the j j j j j=1 j=1 original problem by (31b). c m Define the distance between and given D as d = X j j h ( [k]): j j min j j, and let d = max d . It follows by ^ j j j 2D j j j=1 construction of the update (14) that this distance is halved with each update. To achieve feasibility the distance must be d Thus, p [k] h ( [k]) < 0 for p [k] 6= 0. j j j=1 j made smaller than ( ), which from an initial value of d The last inequality, the line-search condition (23), and the takes dlog e updates. ( ) fact that and are positive, combine to give m m m X X X C. Proof of Lemma 3.1 h ( [k] + p [k]) = h ( [k + 1]) < h ( [k]); j j j j j j Note that j=1 j=1 j=1 0 1 provided p [k] 6= 0. This is equivalent to @ A C x (t) + C x (t ) i 0 i j ` m (0) < m (0) (34) k+1 k j=1 @ when p [k] 6= 0. Consider first that p [k] = 0 for some k 6= 0. = C exp(A(t ))E v ( )d i ` ` ` Then the termination condition (24) is satisfied and statement two holds by (34) and the fact k > 0. Otherwise, p [k] 6= 0 = C E v (t ) i ` ` ` Z for all k 2 N[f0g by hypothesis. In this case, the sequence (m (0)) is monotonically decreasing, and it is bounded + C exp(A(t ))E v ( )d k k=0 i ` ` ` since h are continuously differentiable and each delay must = C E v (t ) i ` ` ` lie in a closed bounded set D . As a result, this sequence converges to a value denoted m i.e., lim m (0) = m < k!1 k C A exp(A(t ))E v ( )d i ` ` ` m (0). This implies that for any > 0 there exists a positive integer K such that m (0) m (0) < ; k K: (35) k 1 k Proof follows from piecewise continuity of v (t ). step d Let < p [k 1] h ( [k 1]). Then j j [k] j=1 d D. Proof of Theorem 3.2 combining the line-search condition (23) and (35) leads to To facilitate the development, relevant variables are indexed < ; k K: by iteration k in what follows. The initial iteration state is [k] k = 0. Each iteration of the SQPA increments the index. Consequently, kp [k 1]k < ; k K; askp [k 1]k 1 1 i) Assume that ( [k]) satisfies (13) and [k] 2 D for j j j j=1 [k]. Hence, the SQPA terminates after K 1 steps and all j 2 N . The updated schedule at k + 1 is [1;m] m (0) < m (0). K 0 [k + 1] = [k] + p [k]; 8j 2 N : j j j [1;m]
Journal
Mathematics – arXiv (Cornell University)
Published: Dec 3, 2020