Two-Scale Model Predictive Control for Resource Optimization Problems with Switched Decisions

Model predictive control (MPC) is widely used in resource optimization problems because it naturally deals with bounded controls and states and allows predictive information to be included. However, at each sampling instant, an optimization problem must be solved. Resource optimization problems with switching control actions naturally lead to optimization problems with integer decision variables, which are computationally costly, particularly when the number of variables is large. As a result, the approach of directly discretizing (DD) the problem to derive a mixed-integer linear program (MILP) sets fundamental limitations on the MPC sampling rate owing to the computational time required to solve the optimization problem. In this paper, we propose a two-scale optimization algorithm (TSOA) for MPC. On the first-scale, the entire prediction horizon is considered and the algorithm provides the optimal resources to be used at each interval with a constant weighting cost. This optimization problem may be cast as a linear program (LP); thus, it is computationally tractable even for a large number of variables and constraints. In the second-scale, the switching nature of the decision variable is recovered by posing an MILP to deploy the optimal resources computed in the previous scale. In this manner, the MILP is solved for a shorter time interval than the entire prediction horizon, thus reducing the number of variables in the optimization problem. The simulation results demonstrate the computational advantages of the proposed algorithm compared to direct problem discretization and optimization.


I. INTRODUCTION
Resource optimization applications with switching decisions are prevalent in electric power systems, such as minimization of the electricity cost by means of load management [1], [2], optimization of hybrid system operation [3], reactive power and voltage control [4], and optimal energy management in liberalized markets [5].
The presence of discrete control variables transforms standard continuous optimization problems [6] into integer optimization problems that are NP-hard [7]; thus, the computational burden increases exponentially with the number of decision variables. Although it is possible to transform optimization problems into mixed-integer linear programs (MILPs) and apply state-of-the-art solvers for planning purposes, the optimization of MILPs becomes infeasible for model predictive control (MPC) applications [8], where the available computation time is limited by the required sampling time. Consequently, the MILP [9] cannot provide the optimal solution before a new sample is acquired [10].
A typical approach to MPC optimization of resource optimization problems with switched decisions is to fix the prediction horizon to 24 h and derive an MILP using a fixed sampling time of one hour, as can be seen from the literature review [1] [2] [3] [11] [12]. This is known as the direct discretization (DD) approach. The DD approach establishes a trade-off between the complexity of the optimization problem to be solved and the suboptimality of the solution owing to large sampling times (e.g. 1 h); however, some applications require incorporating longer prediction horizons beyond 24 h [3], which leads to an infeasible optimization owing to the increase in the problem complexity. Furthermore, even if a feasible solution is found, it has been demonstrated that it may become infeasible when it is periodically repeated [13], thus requiring larger prediction horizons to guarantee feasibility.
Another problem that limits the DD approach in MPC is that it is difficult to determine the computational complexity of the optimization problem to be solved a priori. In fact, as shown in the case study in Section VI, the computational complexity is dependent on the initial state of the system. In this way, the optimization algorithm can provide fast solutions for certain initial conditions, whereas others may require the algorithm to explore a larger number of nodes, dramatically increasing the computational time and rendering the MPC optimization infeasible.
The main line of research to solve the computational feasibility of MPC is to avoid the integer nature of decision variables either by parameterizing the switching times [14] or by relaxations of the optimization problems [14]. In the first case, instead of using binary variables attached to every sampling time, a continuous time description of the switching time instants is used. This approach was successful applied to the control of mixed continuous batch processes [15], control of supermarket display cases [16], and control of wastewater treatment plants [17]. However, the a priori definition of the switching pattern is not known in general and can lead to suboptimal solutions if the pattern is not equal to the optimal one. In the second case, relaxations and reformulations were obtained by adding further constraints and penalization strategies. Unfortunately, although these approaches avoid binary variables, they may lead to nonconvex optimization problems. Consequently, numerical solutions are not guaranteed to be globally optimal.
A different line of research is to tackle the difficulty of the optimization problem using heuristic optimization approaches, such as genetic algorithms, simulated annealing, and particle swarm optimization [18] [19]. These approaches search for the optimal solution using different heuristics. In general, heuristic optimization methods provide suboptimal solutions; however, it is very difficult to establish a priori the required computation time to find these solutions what limits their applicability in MPC.
In this article we present a two-scale optimization algorithm (TSOA) for solving resource optimization problems with switched decisions using MPC. It is derived based on general 1-norm minimization problem that models distinct resource optimization problems under the same formal framework. The major technical difference from previous bibliographical contributions is that TSOA divides the optimization problem into two-scale. In the first-scale, the optimization problem is modelled as a linear program (LP) by integrating binary variables. The LP solution provides optimal resource consumption per interval with a constant cost. In the second stage, switching control action is recovered by solving a MILP for each one of the intervals defined in the previous LP. Second-stage MILPs lead to problems with a smaller number of integer variables than the unique MILP resulting from DD. Furthermore, although the MPC algorithm computes the optimal control action for the entire prediction horizon, only the first control action is applied. As a result, it is not necessary to solve all MILPs for all the intervals defined by the first-scale LP, but only the first MILP that computes the control action that will be applied by the MPC.
The technical contribution of the TSOA is two-fold. First, it is shown that the switched decision problem may be transformed into an LP by integration on an appropriate time scale for the entire prediction horizon, which results in an amenable optimization problem. Second, we define a smaller MILP that recovers the switched control action. Furthermore, the second MILP is designed in such a way that its solution estimates the suboptimality level of the computed solution in a straightforward manner, as shown in Section V. Note that heuristic optimization approaches provide no clues on the level of suboptimality achieved by their solutions.
Summing up, the proposed algorithm can tackle problems with long-time prediction horizons while keeping the complexity of the optimization problem bounded. In this way, the proposed TSOA is scalable with respect to the length of the prediction horizon, because although longer prediction horizons may increase the LP size, the second MILP complexity is kept constant. Furthermore, the TSOA algorithm is scalable with respect to the sampling time selected for discretizing the second MILP because the increase in the number of variables impacts only one MILP (the one that must be solved to provide the first control action for MPC). The results of the TSOA, compared to DD are analyzed in the case study in Section VI. It is shown that the proposed algorithm is a competitive approach for solving MPC problems with switching decisions.

II. FORMAL FRAMEWORK AND PROBLEM STATEMENT
Consider the 1-norm minimization problem [20] as the general modelling framework for resource optimization problems with switched decisions, defined as where norm || · || 1,1 in cost function (1) is a temporal 1-norm and spatial 1-norm, that for a time interval T it is defined as where ||v(τ )|| 1 is the spatial 1-norm of vector v(t) : W x (t) ∈ R n×n and W u (t) ∈ R p×p are the costs or weights in the cost function (1), and are assumed to be diagonal positive semi-definite time-dependent matrices. Constraint (2) is a continuous-time dynamic model represented in statespace form with state vector x(t) ∈ R n , input vector u(t) ∈ R p , disturbance vector d(t) ∈ R q , and dynamic matrix A c ∈ R n×n , input matrix B c ∈ R n×p , and disturbance matrix D c ∈ R n×q . Constraint (4) bounds the state x(t) and constraint (5) sets binary values for control actions. The optimization problem set by cost function (1) and constrains (2)-(5) capture a wide variety of energy and energy cost optimization problems, such as the energy management of a colliery [2], peak load management in steel plants [21], demand-side management of a mine winder system [22], optimization of wind-hydro plant operation [3], and energy cost minimization on pumping stations [23]. Furthermore, although constraints (2)-(5) account for the dynamical system evolution, state constraints, and binary input, it is also possible to consider linear constraints such as production constrains, storage constraints, process flow constraints, sequential constraints, maximum demand constraints, equipment unavailability under maintenance, periodic invariance constraints, and switching costs, among others [13], [21]. It follows that the proposed modelling framework is well suited for the formulation of resource optimization problems. Assumption 1: The state space vector x(t) is positive for all t. Equivalently, the minimum bound in constraint (4) is nonnegative, that is, x min ≥ 0.
Assumption 1 is considered for simplicity only. Otherwise, the LP and MILP optimization problems include new slack decision variables (see [20, p. 294]), which would make TSOA derivation difficult with technical details.

III. OPTIMIZATION PROBLEM FORMULATION IN TWO-SCALE
Instead of solving a single MILP obtained by DD of the continuous-time optimization problem (1)-(5), we define a TSOA approach in which, in the long scale that considers the entire prediction horizon, an LP is solved, and in the short scale, a smaller MILP is proposed that recovers the switching nature of the control actions. The computational benefit lies in the size of the second short-scale MILP, which is considerably reduced from the MILP derived from DD.

A. FIRST-SCALE LP (1 ST LP) 1) LP Time Scale
Weighting matrices W x (t) and W u (t) in the cost function (1) consider the cost of state x(t) and the cost of the control action u(t). For instance, in energy-cost minimization problems, the weighting matrix W u (t) may be the product of the electricity tariff by the electrical power of the actuator. We consider the following assumption for weighting matrices: Assumption 2: Weighting matrices W x (t) ∈ R n×n and W u (t) ∈ R p×p are piecewise constant with K x and K u intervals, respectively. In this manner, the time-dependent weight matrices can be defined as a set of constant matrices W xm for m = 1, . . . , K x and W xn for n = 1, . . . , K u as follows with each W xm and W un being constant diagonal matrices where w xmi ∈ R for i = 1, 2, . . . , n and w ukj ∈ R for j = 1, 2, . . . , p. We define the time scale of the state weight W x (t) as the set of size K x that contains the time instants at which the weighting matrices change, that is T Similarly, we define the time scale of the input weight as If both time scales are not equal, that is, T x = T u , we define a unique time scale T as the union of both time scales with size K; thus, This can be seen graphically in Fig. 1. The goal is to have a unique time scale T, in which the weighting matrices W x (t) and W u (t) are piecewise constants. This time scale T is used to derive an LP optimization problem for the long prediction term.

2) LP Cost Function
In this section, we derive the cost function of the first-scale LP. We begin our transformation by discretizing the cost function (1) with time scale T. First, recall that by definition of the temporal-spatial 1-norm || · || 1,1 , the cost function (1) is Recall that under time scale T, weighting matrices W x (t) and W u (t) are piecewise constants;, hence, equation (12) may be rewritten with constant W xk and W uk as By definition of spatial || · || 1 norm we have that Finally, as the diagonal entries w xki and w ukj of the matrices are non-negative for all k = 1, 2, . . . , K, we have Variables |u j (t)| for j = 1, 2, . . . , p, in cost function (13), are still binary variables, that is, u j (t) ∈ {0, 1}. However, if the binary decision variable u(t) is integrated, the integral ∆t k |u j (t)|dt is no longer binary but real. Hence, we perform the change of variables X ki := ∆t k |x i (t)|dt ∈ R and U kj := ∆t k |u j (t)|dt ∈ R in cost function (13), yielding that is the cost function of a LP with K(n + p) real decision variables X k ∈ R n , and U k ∈ R p , for k = 1, .., K. Finally, cost function (14) can be expressed in compact form as with

3) LP Constraints
To derive the LP constraints, it is necessary to integrate constraints (2)- (5). First, by integrating system dynamics (2) during the time interval ∆t k we have Performing the change of variable D k = ∆t k d(t)dt, and reordering expression (22), the state value at time instant t k is Disturbance D k is known, whereas X k and U k are decision variables. We define the matrix C c := A c B c ∈ R n×(n+p) , equation (23) can be expressed as for k = 1, . . . , K. Applying state constraint (4) to equation (24) yields for k ∈ [1, . . . , K], which can be written in matrix form as where ⊗ is the Kronecker matrix product and Finally, the variable V k := [X k U k ] T (i.e., the integral of the state space and the control action) is also bounded at each time interval ∆t k . Thus, for k = 1, 2, . . . , K, we have which sets the bounds on the integrated state vector X k and on the integrated control action U k .

B. SECOND-SCALE MILP ( 2 ND MILP)
Once the LP defined by cost function (15) and constraints (27) and (29) is solved, the optimal resources per period X * k , U * k , for k = 1, . . . , K, are known. However, the control action in the original problem is a discrete control variable u(t) ∈ {0, 1} p . The discrete control action is recovered by solving the following discrete-time optimization problem for each period with constant weights, that is for k = 1, 2, . . . , K. The objective of the second optimization problem is to determine the optimal binary control u * that best approximates the cost found by the 1 st LP, that is, which best matches control U * k and trajectory X * k . At each period k we state the optimization problem The optimization problem proposed to recover the switching nature of the control action is discrete-time, with sampling time ∆t l , obtained by dividing the interval time ∆t k in L parts, that is ∆t l := ∆t k /L. The cost function (30) minimizes the resources deviation from the LP solution for each interval ∆t k with a constant state and input weights. Equations (31) and (33) are the discrete-time versions of the continuous system dynamics, where the matrices A ∈ R n×n , B ∈ R n×p and D ∈ R n×q are obtained by zero order hold discretization from their continuous time counterparts A c , B c , and D c , with sampling time ∆t l . Discrete-time vectors have the same dimensions as their continuous counterparts, hence x(l) ∈ R n , d(l) ∈ R q , and u(l) ∈ R p .
Remark 1: Dynamic model (31) is the zero order hold discretization of continuous time model (2). Moreover, during sampling times, control action is constant due to its integer nature, as a result both models (2) and (31) yield the same value at sampling times intants ∆t l [24].

1) MILP Cost Function
In the following, we transform the optimization problem defined in (30)-(34) into an MILP. The cost function (30) can be written in matrix form as follows min u(t) ∆t l W xk I nxn I nxn . . . I nxn . . .
which may be transformed into an MILP by introducing slack variables r ∈ R n×1 and t ∈ R p×1 , see for instance [20, p. 294], yielding

2) MILP Constraints
Now, we transform constraints (31)-(34) to be used in MILP. Applying the dynamical system in (31) to state constraints (33), for l = 1, . . . , L, we obtain the following equations in matrix form: with variables and prediction matrices defined as VOLUME 4, 2016 Finally, the constraints in (38)-(39) can be compacted as

IV. TWO-SCALE MODEL PREDICTIVE CONTROL ALGORITHM
In this section, we describe the two-scale optimization algorithm used in MPC. First, on the long scale provided by the time scale T = {t 1 , t 2 , . . . , t K }, an LP is defined and solved, as shown in Figure 1. Therefore, this approach is feasible for long prediction horizons. As a result, the 1 st LP (summarized in Section IV-A) provides, at each time interval ∆t k with constant weights, the optimal resources V * k = [X * k U * k ] T . Next, for each short time scale defined by the time interval ∆t k , where k = 1, 2, . . . , K, an MILP (summarized in Section IV-B. 2 nd MILP) is solved to determine the optimal switched control u(t l ) * . In Figure 1, we shown the 2 nd MILP for ∆t k and ∆t k+1 . The 2 nd MILP is now feasible because the time scale is reduced to ∆t k , instead of T , as is the case with the MILP obtained by the DD of problem (1)-(5). Furthermore, for the optimization in MPC, only one MILP, the one whose its ∆t k contains the present time instant, is solved to obtain the optimal control action to be applied by MPC.
The1 st LP is defined on the large time scale T = {t0, t1, . . . , t K }. It provides the optimal resources X * k and U * k for each time interval ∆t k . The time intervals ∆t k := t k − t k−1 for k = 1, 2, . . . , K of time scale T define the intervals where the 2 nd MILPs are solved. In particular we shown the 2 nd MILP for ∆t k and ∆t k+1 . However, for the optimization in MPC only one MILP, the one whose ∆t k contains the present time instant, is solved to obtain the optimal control action to be applied by MPC.
Once the fundamental idea behind the two-scale algorithm is formulated, we describe its implementation as an MPC control problem. Figure 3 shows the flowchart of the TSOA MPC implementation. After loading the problem data and the initial value of the state x 0 , the 1 st LP time scale T is computed and the 1 st LP problem formulated. The solution of the 1 st LP is given by V * k = [X * k U * k ] T , that is used to formulate the 2 nd MILP. Its solution provides the vector u(t l ) * for l = 0, 1, . . . , L − 1. Then, the first control action u(t * 0 ) is applied, the new state x 0 is measured, and the process is repeated as shown in the flowchart in Figure 3.
The major advantage of the two-scale optimization algorithm for its application to MPC lies in the fact that it is not necessary to solve all the 2 nd MILPs for ∆t k with k = 1, 2, . . . , K but just one 2 nd MILP for the current ∆t k , that is, for ∆t k that contains the present time instant. Although the two-scale approach provides computational advantages for planning by reducing a single large MILP into K-shorter MILPs, its major advantage relies on its application to MPC because only an LP followed by a single MILP on a short scale must be solved to compute the optimal control action to be applied at each sampling instant.
In fact, for a time horizon T = K∆t k , a sampling time ∆ t l = ∆t k /L, and a dynamical system of order n and p inputs, the DD approach results in a MILP with n(KL + 1) + p(KL + 1) variables and 2nKL constraints. On the contrary, the TSOA results in a MILP with n(L+1)+p(L+1) variables and 2nL constraints. As a result, the DD yields a problem which number of variables and constraints are of order O(KL) whereas in the TSOA are of order O(L). As the MILP is an NP-hard optimization problem, the complexity increase is not linear in K but factorial, thus dramatic reduction in computation time can be achieved even for moderate values of K. In Section VI, we quantitatively demonstrate the computational advantage of the TSOA in MPC when compared to the computational cost of the DD approach.

A. FIRST-SCALE LINEAR PROGRAM
This section summarizes the 1 st LP cost function and constraints. These are

B. SECOND-SCALE MIXED INTEGER LINEAR PROGRAM
This section summarizes the 2 nd MILP cost function and constraints. These are (56)

V. OPTIMALITY AND COMPLEXITY PROPERTIES OF THE TWO-SCALE ALGORITHM
In this section we analyse the optimality and complexity properties of the proposed two-scale optimization. First, we define the optimal solution of problem (1)-(5) as u * * (t) and x * * (t). This baseline serves as a comparison to other solutions. We also define the optimal solution of the 1 st LP (44)-(48) as X * k , U * k for k = 1, 2, . . . , K, and the optimal solution of the 2 nd MILP (49)-(56) as u * (t) and x * (t).
We have the following lemma that bounds the optimal solution u * * (t) and x * * (t) with the solutions of the two-scale algorithm: Lemma 1: The optimal solution of the optimization problem (1)-(5) is bounded above by the optimal solution of the TSOA and below by the optimal solution of the 1 st LP, that is, Proof 1: Consider that the optimal solution of optimization problem (1)-(5) is given by x * * (t), u * * (t). As a result, it follows from the definition of optimality that ∀ x(t), u(t) Particularizing the previous result for x(t), u(t) the solution provided by the TSOA algorithm proofs the first inequality. The second inequality follows because the 1 st LP is a  relaxation of the optimization problem (1)- (5). In fact, the constraint of problem (1)-(5) is The particularization of t ∈ R into time instants t k , k = 1, 2, . . . , K and the integration of x(t) for each ∆ k yields the following constraints of the 1 st LP, VOLUME 4, 2016 As a result, it is not true that equations (60) and (61) implies equation (59).
The main advantage of Lemma 1 is that the TSOA automatically provides the upper and lower bounds on the optimal solution of the problem. As a result, if both bounds are tight, we can ensure near-optimality properties of the solution. In the following, we show that the cost function of the 2 nd MILP automatically provides an upper bound on the suboptimality of the two-scale algorithm.
The optimal solution of the 2 nd MILP at interval ∆t t is x * (t l ) and u * (t l ) for l = 1, 2, . . . , L. Furthermore recall that, for any norm, we have the inequality ||a − b|| ≥ ||a|| − ||b|| . As a result, it follows that The term (62) is the discrete-time approximation of the solution computed by the two-scale optimization algorithm for the problem (1)-(5) because whereas the term (63) is the optimal value obtained by the 1 st LP. As a result, the cost of the 2 nd MILP is greater than or equal to the difference between the cost of problem (1)-(5) obtained by the solution provided by the two-scale algorithm minus the cost of the 1 st LP.
Lemma 2 allows to obtain optimality bounds on the two-scale solution. In fact, Lemma 2 can be rewritten as, where J k is the 2 nd MILP optimal cost, J T SOA k is the optimal cost obtained by the two-scale optimization algorithm, and J LP k is the 1 st LP optimal cost, all considered at time interval ∆t k . As a result it follows that If these bounds are tight we have a certificate of nearoptimality for the TSOA for each interval ∆t k .

VI. CASE STUDY
In this section, we consider the optimal operation of a water supply system as a simple benchmark to show the infeasibility of MPC arising from MILP optimization problems derived by direct discretization from optimization problem (1)-(5). The benchmark comprised three reservoirs and two pumps to be controlled, as shown in Figure 4. Even for this simple problem, infeasible computation results were obtained for the MPC with DD optimization.
First, we demonstrate how infeasibility can be attained by optimization of the DD model. Consider optimal pumps MPC for a 24 hours time horizon with initial state x 0 = [200 100 100] T m 3 . The optimal problem is transformed into an MILP with a sampling time T s = 30 min. The solution was computed using the 'intlinprog' function in MATLAB with a computational time of 0.4364 s. However, when the initial state is changed to x 0 = [100 30 30] T m 3 , the optimal solution is not found, because the 'intlinprog' function reached the default value of the maximum number of nodes with a total computational time of 1815.04 s (i.e.

min).
A consequence of the previous example is that it is very difficult to establish a priori the computational requirements to find a (sub)optimal solution with an MILP obtained using DD. In fact, it is the change in the initial condition x 0 that transforms an easy optimization problem into a difficult one, by increasing the required computational time by two orders of magnitude. This is because the value of the initial condition, x 0 , defines the problem constraints. Thus, for certain values of x 0 it may be the case that a larger number of constraints are active, increasing the required computational time to find the optimal solution. Furthermore, note that this approach is infeasible for use in MPC because the required computational time of 30.25 min is larger than the sampling time of T s = 30 min. As a result, the direct approach of using DD in MPC is, generally, not a feasible approach, even for problems with a reduced number of actuators and constraints, and with large sampling times. The fact that one instance of the optimization problem can be solved within the allowed computational time, does not guarantee that the next instance of the optimization problem will be solved. Now, we perform a quantitative comparison analysis between DD approach and the TSOA for the initial state x 0 = [200 100 100] T m 3 . First, we obtained the optimal solution by DD with two distinct sampling times of 30 min (i.e., DD 30 ) and 5 min (i.e., DD 5 ). The cost of the direct discretization algorithm with a sampling time of 30 min is DD 30min = 195.85 whereas the cost with sampling time of 5 min was DD 5min = 166. 22. The reduction in the sampling time from 30 to 5 min reduces the cost of 17 %; however, the required computational time was increased to 93 %. This clearly shows the trade-off between optimality and computational complexity in the DD approach.
For comparative purposes, we now compare the DD approach with the TSOA, where the 2 nd MILP is discretized using the same time interval as in the DD approach. The optimal cost achieved by the TSOA approach is equal to that obtained by DD, as shown in Table 1. The total computational time (i.e. 1 st LP time plus the 2 nd MILP) for TSOA 30 is a 23 % larger than the DD 30 approach. However, for the TSOA 5 the computational time is 25 % shorter than the DD 5 . Furthermore, the TSOA computational time hardly increases with a reduction in the sampling time from 30 min to 5 min, showing the scalability properties with respect to the sampling time of the proposed two-scale optimization approach.
The TSOA also provides optimality bounds for the computed solution. First, recall that from Lemma 2, we have J LP k − J k ≤ J T SOA k ≤ J LP k + J k . In this way, for the first time interval ∆t 1 and sampling time 30 min, we have J LP1 =138.48 and J 1 =21.76; thus, the cost of TSOA 30 is bounded by 116.72≤ J T SOA1 ≤160.24. Now, applying Lemma 1, the optimal solution for the whole prediction horizon is bounded above by J T SOA =195.85 and below by J LP =158.26. As a result, the a priori worst-case suboptimality of the TSOA with respect to the DD is approximately 23 %. A posteriori, it can be seen that both the TSOA and DD, provide the same cost. For a sampling time of 5 min we have J T SOA =162.22 and J LP =158.26. Consequently, the a priori worst-case suboptimality of the TSOA with respect to DD is 2.5 %. Henceforth, the provided TSOA optimal solution, in the worst-case, deviated by less than 2.5 % of the optimal solution. Now, consider the more computationally demanding MPC problem with the initial state x 0 = [100 30 30] T m 3 . The DD 30 cost is 515.79 whereas the DD 5 cost is 460.35.  Reducing the sampling time from 30 min to 5 min reduced the cost by 10 %. However, although DD 5 has a computational cost of 1.08 s, for DD 30 , the computational cost is 1855 s because, although DD 5 has more variables than DD 30 , the shorter sampling interval makes it easier to find a feasible solution, whereas DD 30 has to explore more nodes to find a suboptimal solution. This again shows the difficulty to ascertain the computational complexity of the DD problems. We now compare the TSOA approach with the DD. TSOA 30 achieves a cost 2.3% lower than DD 30 , as TSOA 5 and DD 5 achieve nearly the same cost, as shown in Table 2. Regarding the total computational time for TSOA 30 it is 0.38 s and for TSOA 5 it is 0.39 s. The computational time is shorter than that in the DD approach, and it can be seen that the computational time hardly changes with respect to the sampling time, again showing the scalability properties with respect to sampling time.
Regarding the optimality bounds for the first time interval ∆t 1 and sampling time 30 min, we find that J LP1 =276.96 and J 1 =9.89; thus, the cost of TSOA 30 is bounded by 267.07≤ T SOA 1 ≤286.85. For the first time interval ∆t 1 but with a sampling time of 5 s, we have J LP1 =276.96 and J 1 =0; thus, the cost of TSOA 5 is bounded by T SOA 1 =276.96, that is the cost achieved by the TSOA 5 . Applying Lemma 1, the optimal solution for the whole prediction horizon is bounded above by J T SOA =503.47 and below by J LP =457.63. As a result, the a priori worst-case suboptimality of the TSOA with respect to DD is approximately 10 %. A posteriori can be seen that both, TSOA and DD, provide the same cost. For a sampling time of 5 min we have J T SOA =459.23 and J LP =457.63. Consequently, the a priori worst-case suboptimality of the TSOA with respect to the DD is 0.3 %. Henceforth, the provided TSOA optimal solution, in the worst-case, deviated by less than 0.3 % of the optimal solution.
Finally, we graphically show the optimal control and state evolution for the initial state x 0 = [100, 30, 30] T , solved with DD 5 and T SOA 5 . The optimal control actions for DD 5 are shown in Figure 5, and for T SOA 5 in Figure 7. The state evolution, that is, the reservoir levels, is shown for DD 5 and T SOA 5 in Figure 6 and 8, respectively.

VII. CONCLUSION
In this article, a new two-scale optimization algorithm (TSOA)is proposed for use in model predictive control (MPC) for resource optimization problems with switched decisions. The proposed approach overcomes the infeasibility and scalability problems in optimization that appear in the mixed-integer linear program (MILP) derived from the direct discretization (DD) approach because the TSOA can consider the entire prediction horizon while keeping the MILP complexity bounded. As a result, the solution is feasible as the computational time is reduced with respect to DD, and is scalable because the increase in the prediction horizon does not increase the MILP complexity. Furthermore, the algorithm is also scalable with respect to the sampling time because only one MILP with a fixed length must be solved to compute the first control action. The TSOA was tested on a three reservoirs with two pumps water supply system. The computational results show that MPC with TSOA is computationally feasible because the required computational time is shorter than that of the DD approach. Furthermore, it is scalable as the computational time hardly increases when the MPC sampling time is reduced. As a result, the TSOA enables the application of MPC in optimal resource problems with switched decisions.