Hierarchically coordinated economic MPC plantwide control of mixed continuous‐batch units in process industries with application to a beet sugar plant

This paper deals with the optimal operation of processes that combine batch and continuous units. The proposed approach is based on a hierarchical architecture with two layers. The lower layer manages the optimal operation of the batch units in a noncentralized fashion, whereas the upper layer is concerned with the scheduling of the batch units and with their smooth integration with the continuous ones. The strategy has been applied to a challenging problem found at the interface of the fed‐batch crystallizers with the upstream evaporation section in beet sugar factories. The presented approach aims at designing a system that can effectively perform the plantwide economic optimization, specifically addressing the issue of energy efficiency, while complying with all process and operational constraints.


INTRODUCTION
Traditionally, the plantwide management of process industries has been undertaken using a decentralized control approach, in which each main unit of the factory's flow sheet was controlled individually by single variable regulation loops of temperature, pressure, level, composition, and other conventionally measured variables. In the best case, some advanced multivariable schemes were put in place to keep the operation in a safe and reasonably efficient working regime. The material and energy couplings between these controlled units were treated as perturbations whose negative impact was to be kept in check by the right control design, including an appropriate pairing of controlled variables with a controller's manipulated outputs and with a conservative process design that included, for example, the placement of intermediate buffer or surge tanks. The plantwide management was more an art, refined by the experience of process and control technologists, than a science backed by the academia, a situation denounced in the work of Foss. 1 There were general rules proposed such as those in other works, 2-4 which remains valid.
The situation has changed for the better in the last decades, a change driven by the pressing need to be competitive in a global economy. The design of process factories moves to a greater heat integration, with an increasing presence of recycling streams and the reduction of inventories leading to a more tightly coupled plant. The plantwide control has moved in accordance. It is not enough to have a smooth, free of bottlenecks, and reasonably well controlled plant.
Ideally, the plantwide control task is to choose the degrees of freedom that might remain after complying with the needed constrains, in a rigorous, economically optimal manner, solving problems such as min. (1) s.t. h plant ( .
In (1)(2)(3), the plantwide operation is cast as a mathematical programming problem. There is an economic performance index (J), eg, the factory's profit, that could be stated as a balance between the cost of producing the marketable product and its market value. The instantaneous profit or stage cost ( e ) should be totalized over the whole considered exercise by performing the integration indicated in (1). The optimization is to be performed choosing the plantwide degrees of freedom represented by u g . In a natural setting, the parameters p cost and p value represent the cost of raw materials and utilities, such as energy on one side and the price of the final products on the other, whereas W cost and W prod are the respective flow rates. The optimization problem is constrained by the function h plant , which stands for the time (t)-dependent dynamic-algebraic (differential algebraic equation [DAE]) model of the plant with dynamic state x, algebraic variables z, and parameters p. The inequalities g plant establish the technological, safety, environmental, and other constraints that might apply. All the previous elements can be considered as column vectors in the general case or vector-valued functions. The ones conforming the scalar index J should have the appropriate dimensions as to render the represented scalar products well defined. The plantwide solution of (1-3) is far for being trivial. The dominant paradigm in which the problem is framed is based on a layered approach (Figure 1) that reflects a natural division of concerns and the different timescales involved. 5,6 It stretches from the economic and long-term planning and scheduling of the factory at the very top of the hierarchy in the enterprise resource planning (ERP) layer as codified, for example in the ISA-95 standard, 7 all the way down to the basic regulatory control loops and then to the process itself at the bottom ( Figure 1A). The link between the ERP business-related level 4, which works in a time frame of days, weeks, and even months, with the supervisory and advanced control level 2 that should be able to respond to the faster dynamics of the underlying process, is the task of the level-3, manufacturing execution system (MES) layer. In Figure 1B, the distribution of responsibilities of a typical plantwide control systems FIGURE 1 ISA-95 layered approach to plantwide management (A) and the control pyramid (B) follows also a layered approach. The planning and scheduling of the control stack is dedicated to provide short-term, more operative actions belonging conceptually to the abovementioned MES layer. On top of the regulatory layer of Figure 1B, the advanced control function is typically implemented by the known model predictive control (MPC) controllers 8,9 of widespread use in process plants. They can handle the complex couplings in multivariable problems. They are provided with a dynamic mathematical description of the plant used to optimize the operations according to a given performance index. The optimization setup can easily incorporate operational constraints and inputs limits, being this one of its most attractive features. The mathematical optimization is carried out in an open-loop fashion for a finite number of sampling instants into the future. Disturbances are taken into account by means of the celebrated receding horizon strategy. In this context, MPC is of the tracking or regulatory type with a weighted quadratic index that penalizes the deviation from the set points and the control effort. The manipulated variables of this layer are handed down as references to be followed in the regulatory layer while the set points to track are provided by the real-time optimization (RTO) layer directly above. The RTO 10 solves another optimization problem, this time with an economic objective and a static, accurate, typically nonlinear model of the factory. To do that, it must wait for the plant to be at steady state, a clear drawback that limits its applicability. Other criticisms of the described scheme point out the difficulties of having to keep up to date two different models, with different characteristics, at the RTO and MPC layers, which however need to be consistent enough so that the demanded economical best regime by the upper layer remains feasible and implementable in the lower one. With legitimate criticisms apart, the pyramidal structure just described, suitable to be embedded in hierarchical structures as the one proposed by standards like ISA-95, constitutes a useful paradigm.
An alternative is represented by the so-called economic MPC (EMPC), where the quadratic set point-following objective function of the tracking or regulating MPC is substituted by a free-form index, which is typically related to the economic performance of the plant. This apparently simple change, which in practice conflates into a single problem the RTO and MPC layers above, has profound consequences. On the plus side, the problem of the existence of two models is circumvented as is also the need for waiting for the steady state of the plant. The promise of EMPC is to make the transitions along the unavoidable dynamic evolution of the plant in an economically optimal manner. The challenge resides in achieving the stability and performance requirements and in the nonlinear nature of the problem. The considerable literature dealing with stability for the conventional MPC case, as discussed in the survey, 11 were based of the special form of the index function. The new free-form index of EMPC calls for another type of treatment as discussed, for example, in the tutorial 12 and treated in detail in an already important number of contributions. [13][14][15][16][17][18] However, a monolithic implementation of EMPC hardly constitutes a valid alternative for plantwide, large problems. The difficulties here not only reside in the prohibitive burden of computationally finding the solution every sampling time but also on the rigidity and maintenance difficulties of a single, all-encompassing mathematical model. Simply put that this monolithic approach is not the way that control solutions for plantwide problems are developed and grow over time.
Therefore, closing the circle, there is a growing trend to develop, not the completely decentralized schemes of the past, but distributed solutions such that, upon partitioning the plantwide problem in some suitable way, develops controllers to each specific subplant while maintaining a good performance for the whole plant. In general, the fact that the different subplants are not independent of each other demands the interchange of information between controllers, either in a peer-to-peer basis, or using an upper layer coordinator to take care of the connecting physical variables.
In this contribution, we put forward some ideas to tackle a loosely coupled, hierarchically coordinated EMPC solution as applied to a problem where there is a group of parallel units sharing a set of common resources, combining continuous and batch processes. The novelty of the approach is that the competing units are of a batch or fed-batch type, so the optimal operation has to embed into it the solution of the scheduling problem for the discontinuous units, coordinating their operation with the continuous subplant.
The strategy developed is applied to a case study representing the interface between the continuous upstream department of a beet sugar factory and the fed-batch crystallizers. The optimal problem is cast in a way that the plantwide solution hinges on the concept of sugar factory energy efficiency.
The simulated case study used is described in the work of Mazaeda et al 19 and the solution adopted extends the previous work reported in conference proceedings. 20,21 This paper is organized as follows. Section 2 goes deeper into the problem of distributed optimal control, to establish, in a general framework, the characteristics of the type of problem specifically addressed here and the solution adopted. Section 3 describes the beet sugar plantwide case study, highlighting the main tradeoffs to be faced in its optimal management. Section 4 cast the solution to the sugar plant in the framework previously exposed. Section 5 discusses the obtained results and their relevance in terms of energy efficiency.

COORDINATED OPTIMAL CONTROL OF FED-BATCH PARALLEL UNITS SHARING COMMON RESOURCES
There is a large literature in distributed peer-to-peer MPC, classifying the different alternatives according to how the information is interchanged among the different control agents, the nature of the function to optimize, the nature of the variables of the underlying process that are shared, the linearity or not of the model used, if iterations are required or not at every sampling instant, among other factors. For example, it is useful to classify noncentralized solutions as cooperative if the control agents share a common performance index and noncooperative otherwise. Relatively recent accounts are given in the works of Christofides et al 22 and Scattolini. 23 Some algorithms exhibit proof of stability and/or convergence to the plant optimum, typically under assumptions of linearity or convexity of the model and constraints. An important amount of work has been dedicated to the tracking distributed MPC (DMPC). Recent contributions, for example, deal with hybrid models 24 and event-triggered schemes. 25 In the works of Maestre et al 26 and Alvarado et al, 27 different algorithms are compared as applied to nontrivial benchmarks.
Distributed EMPC (DEMPC) is also an active field of research; see for example other works. [28][29][30] The use of individual controller agents to solve the plantwide problem is predicated upon the existence of a reasonable partition of the original model. There are contributions dedicated to the automatic division of the plantwide problem. 31 The objective is to obtain a number of clusters, which are internally highly coupled and with scarce connection among them. In many cases, however, the topology of the process factory, where units of a certain complexity and a very well defined functionality are connected by product material streams and energy interchange paths, suggest a natural distribution of loosely coupled tasks.

Hierarchically coordinated optimization: Problem statement
An alternative to the peer-to-peer distributed control with interchange of information between the controllers is the use of a hierarchically coordinated approach. The coordination is carried out by a lightweight upper layer, which is in charge of receiving information from the subordinate controllers and of adjusting them with the purpose of achieving the desired optimal plantwide behavior.
The use of divide-and-conquer strategies to optimize large problems is not a new subject. 32,33 An influential article in applying this kind of scheme in process industries is that of Jose and Ungar. 34 It advocates an optimal distributed solution based on dual decomposition or price coordination but cast in a certain way by means of carefully designed slack variables. This kind of solution constitutes a type of Lagrangian relaxation stemming from the known literature on optimization. 35 An algorithm based on price auctions would arrive iteratively at the solution of the overall problem. Although Jose and Ungar 34 give a microeconomic interpretation of the problem in the form of a supply chain metaphor applied to the relation between the upstream and downstream process units; it is worthwhile to point out that the prices used in this context are simple devices to reach the consensus among the controllers. Another application dealing with price coordination of oxygen networks, this time stated in the EMPC framework, which used a PID inspired procedure for adjusting prices, is given in the work of Martí et al. 36 The procedure implies, as is often the case, the use of iterations on every sampling time to arrive to the consensus solution.
The specific contribution offered here is applicable to the similar case where different units in parallel compete for a set of shared scarce resources. This type of topology is not an exception because, as commented in the works of Downs and Skogestad 3 and Jäschke and Skogestad, 37 it is the pragmatic alternative for increasing throughput of an existing process plant. The complicating characteristic and the novelty of the approach are that all or most of the units are of batch or semibatch character. The proposed hierarchically coordinated EMPC approach should steer the plant to an optimal or at least a better overall economic regime solving the embedded scheduling problem introduced by the presence of the discontinuous units. Due to this fact, in this continuous-batch scenario, the execution of iterations at every sampling instant, in a situation where the shared resources are tightly constrained, is an untenable proposition. The seeking for convergence repetitions would need to be reiterated, on its turn, for the whole horizon into the future for each of the parallel batch units and repeated all over again, possible many times, until arriving at the optimal scheduling. Therefore, we assume in this case that the common resources can fluctuate between limits without putting in jeopardy the independence of the subproblems sharing them. In the process industry interpretation, this means that we assume the existence of buffer containers where the inventories of the shared resources, expressed for example as the levels attained by liquids in tanks or by the pressure of gasses contained in closed vessels, are free to move between certain limits without compromising the independence of the units using them. Note that this is not an unreasonable assumption because this is the almost universal practice in this kind of situation. In distributed optimal literature when the subtasks do not have to agree on the values of some common connecting or troublesome variables, the problem is classified as embarrassingly parallelizable. It can be said that here we pragmatically choose to keep the problem that way by enforcing some added restrictions on the common inventories. In any case, an optimal solution to this kind of problem with the right scheduling included, in addition to all other operational and economic benefits, would allow for a reduction in size of the intervening buffer tanks.
Therefore, we set out to solve the optimal problem stated in (1-3) for the kind of hybrid plant with parallel batch units above described. The problem is posed as a two-level hierarchically coordinated optimization. The plantwide problem is considered to be implemented, in the general case, with EMPC controllers for the coordinator and for the subordinates. The upper-layer coordinator interchanges information with each one of the controllers of the subordinate units. The coordinator problem (4-7) is provided with the plant real economic index (4), which is identical to the original one (1). The equality vector-valued function h coord in (5) defines a very simple model basically describing the evolution of intervening inventories of the shared resources. The vector of inequalities g coord (6) makes sure that the limits of the inventories that keep the subproblems noninteracting are satisfied. The decision variables include the vector of starting commands (Load) for each one of the subordinate batches for the whole considered horizon into the future. Similarly, it could optimize over some other degrees of freedom here no specified and represented by u g . The vector constitutes a device for performing the coordination with the lower-layer controllers and can be interpreted as normalized prices. They are constructed as real values in the closed interval [0,1] and represent, as going to be clarified later on, the relative importance that each particular subproblem gives to the rate of production of marketable stream over the cost of the consumption of utilities as for example energy. The variables are the knobs that the coordinator has for steering each subordinate controller to achieve the objectives of the plantwide problem. However, to do that, the coordinator needs information from the lower-layer controllers: It needs the current, updated profile of consumption of the shared resources (w shared ) and the sensitivity of those profiles to the mentioned tuning knob variables over the whole batch. Note that this profiles are, in general, time varying due to the semibatch nature of the subordinate processes. That information is needed to construct a linear approximation of the effect of the mentioned artificial prices on the shared inventories. For the linear approximation to be valid, a restriction is added (7), where vector prev are the prices used in a previous iteration and the α parameter vector, with small-enough elements, is used to define an appropriately narrow trust region. The operation indicated in (7) is the Hadamard product of the corresponding vectors.
Each one of the batch-type subordinate optimal problems is defined as in (8)(9)(10). Note that the problems are posed in terms of their own cycle time ( i ), anchored at the beginning of each batch and which is displaced with respect to the real time (t) of the coordinating problem according to the implemented schedule. The decision variables (the vector u local ) are privative of each local problem, and the same happens with the dynamic and algebraic variables represented by the vectors x local and z local . The shared variables (vector w shared ) are independent if the restrictions (6) of the coordinator layer are met. The h local and g local vector-valued functions represent the DAE equations defining the subordinate unit's dynamic evolution and the path restrictions that need to be enforced respectively. min.
The local decision functional to optimize (8) leads to a maximization of a kind of profit defined at the local level. It encourages a higher production of the intended product W prod with the smaller possible consumption of the needed utilities W cost at the shortest possible batch cycle time (T batch ). In addition, this is what would happen if the real market prices of product and cost were used in (8). However, this need not be the best overall behavior: Nothing is gained by conducting each parallel batch units at a higher-than-needed pace because the overall economic performance is the one that counts. It is true that the internal logic of an economically sound process factory design usually implies an optimal operation that maximizes throughput, until the plantwide bottleneck is hit, as discussed in the work of Acebes et al. 38 However, a particular batch unit might or might not be in that bottleneck, and that statement is true even for the combined effect of all the parallel batch processes because the factory rate might be fixed by some other upstream or downstream unit. Note that in the case where the production rate is decided elsewhere, ie, by the upstream plant, for example, in the case of a partial plant, the economical best line of action is to keep up, in average, with the already decided product stream flow rate, maximizing the efficiency by reducing costs, eg, the consumption of energy.
In the light of the previous discussion, the batch subordinate problems use, instead of the real prices, the already mentioned normalized parameter that is going to be adjusted incrementally by the upper-level coordinator, that, from its vantage point, has a whole picture of the plant economics. Therefore, for the sake of the classification of the proposed hierarchically coordinated approach, it can be said that it belongs more to the class of cooperative algorithms 23 although the original definition does not apply strictly because the function to optimize is not exactly the same. Note also that the coordination performance function needs not be the sum of the ones corresponding to all subordinate units. The problem could best be classified as belonging to the hierarchically coordinated noncentralized problem with consistency of interest as defined in the work of Findeisen. 33

Hierarchically coordinated optimization: Solution by means of EMPC
The coordinated problem of the previous section can be solved in the framework of hierarchically coordinated EMPCs ( Figure 2): one for solving the master problem and as many subordinated controllers as there are parallel units ( Figure 2A). The subordinate semibatch EMPCs are implemented in the usual shrinking horizon scheme using a sampling frequency suitable to each individual problem dynamics. The subordinate EMPC should be ready to respond to the requests of the coordinator by giving to it the most recent complete predicted profile for the consumption of shared resources along with the sensitivity of those profiles to the local parameters. The coordinating EMPC has in general a different sampling frequency, and must construct on the fly, with the information obtained from the subordinate controllers and other measurements, and basic mass and energy balance equations, a simple model of the relevant inventories. The prediction and control horizons for the coordinator are chosen as an integer number of complete cycles into the future for each of the subordinate units. The scheduling problem would suggest the need to solve, at each iteration of the coordinator, a mixed-integer nonlinear problem (MINLP) that could be very demanding in the computational sense. In this particular case, this could be avoided by means of the parametrization described in the work of de Prada et al, 39 where the Load i binary command for starting subordinate cycles is replaced by the continuous time T loadi to wait from the considered time to the instant of batch initiation as shown in Figure 2B.

PROCESS DESCRIPTION
The hierarchical approach described in the previous section would be applied to a realistic sugar industry plantwide management problem. To best understand it, a digression for describing the chosen case study follows.
Beet sugar factories are conventionally divided in two main sections: the beet and the sugar departments. The upstream beet section in charge of extracting the sucrose from the sliced beets in a diffusion in water process, of eliminating the best part of the impurities that were extracted along with the sucrose in a purification section and, finally, of concentrating the resulting solution by means of a series of evaporators. The beet department, which is an essentially continuous plant, delivers the resulting syrup, a technical solution of sucrose plus impurities in water, to the downstream sugar department, which is in charge of producing the commercial sugar grains. The sugar house consists, typically, of three stages of similar structure. The first stage produces the commercial sugar grains. It consists of an array of parallel crystallizers in the form of fed-batch units. They consume the feeding liquor provisionally kept in a buffer tank to deliver, at the end of the cycle, the so-called massecuite, ie, a slurry where the fully developed sugar grains are embedded into the mother liquor, ie, a solution containing almost all the impurities fed to the unit along the batch and the dissolved sucrose that has not migrated to the crystal phase. The massecuite is provisionally kept in tanks in its way to the centrifuges where the sugar crystals are separated from the mother liquor. The mother liquor of the first stage has still too much sucrose for it to be thrown away, so it is further processed in the a second and then a third stage. All the sucrose recovered, in the form of sugar crystals, from the second and third stages and other intermediate syrup streams are returned to the feed liquor tank of the first stage, a strategy that is essential to secure the economic viability of the plant.

Plantwide case study description
The chosen case study considers the last three effects of a multiple effect evaporation section and three identical fed-batch crystallizers placed in parallel ( Figure 3). A key element is the floating level buffer tank used for accommodating the continuous supply of concentrated syrup from the last effect of evaporation with the intermittent and time-varying combined consumption of the crystallizers.
The evaporators and the crystallizers are tightly connected. Evaporators not only deliver the stream of the syrup to process but they also constitute the immediate source of energy in the form of steam which is required by the crystallizers. The pattern of syrup and steam intake of each crystallizer present a time varying profile along the cycle, a behavior which is typical of fed-batch equipment. The crystallizers compete for the common resources represented by the availability of steam and of syrup. Therefore, the plantwide problem comes with the need of scheduling crystallizer's cycles in an appropriate manner.
To better understand the problem, a brief discussion of the workings of the evaporation station and of a fed-batch sugar evaporative crystallizers are in order.

Evaporation station
The evaporation section is organized in a multieffect scheme that allows for an efficient reutilization of the primary source of steam produced at the factory's boilers. Each individual effect works in a continuous fashion, evaporating a fraction of the water of its incoming juice. The energy to achieve that is given by the latent heat of condensation of the steam fed at the primary of an incorporated heat exchanger. The heating steam condenses in the exterior shell and communicates heat energy to the juice at the other side of heat transfer wall. Part of the water evaporates from the boiling juice, and ascends to the top of the chamber, whereas the remaining, more concentrated juice goes out of the level-controlled effect. The diluted thin juice coming from purification enters the first effect, which also receives the live steam from the boilers. However, each successive unit downstream of the first, takes the more concentrated juice from the previous one but, and this is key for the reutilization of energy, it uses as energy source also the vapors evaporated from its upstream unit. For this scheme to work, the boiling process in each effect has to be conducted at a lower pressure, and thus temperature, as compared with its upstream neighbor. The final effect delivers the concentrated juice or syrup while its boiling pressure, typically below atmospheric pressure, is fixed by a barometric condenser. The degree of concentration (or brix) of the delivered syrup is regulated by adjusting the flow of vapor to the condenser and thus, and in conjunction with the live pressure value at the input of the station, it determines the boiling pressures of all the intervening effects. Not all the vapor produced at each effect is used in providing energy to the downstream one. A considerable part of the flow is diverted to supply energy to different units of the factory. One of the most important consumers of the reused energy that the evaporation delivers are the downstream crystallizers.

Sugar fed-batch evaporative crystallizer
The management of a sugar industrial crystallizer must pay close attention to the supersaturation of the solution (Figure 4). Supersaturation is defined as the ratio between the current concentration of sucrose in the solution and the one that defines the solubility at the current temperature and purity. The supersaturation should be kept at moderate values along the whole process, ie, safe into the so-called metastable zone ( Figure 4D). An infrasaturated solution leads to the dissolution of existing crystals, whereas a too-high supersaturation risks entering the labile zone where spontaneous uncontrolled nucleation of new crystals might affect the crystal size distribution (CSD) of the population of sugar grains, rendering the batch useless in the more serious cases. Supersaturation can be created in a number of ways, but in the case at hand, it is done by evaporating the required amount of water.
The sugar fed-batch crystallization process is conducted inside a vessel where the needed evaporation of water is carried out at vacuum conditions to keep the boiling temperature low enough as a way to avoid the thermal degradation of sugar. The unit ( Figure 4A) has an incorporated heat exchanger, similar to the ones described above for the evaporators. Both the steam pressure at the primary of the heat exchanger and the vacuum inside the chamber are regulated by their respective PID loops. The process is conducted by following a sequence of steps or recipe. As shown in Figure 4B, the cycle starts by introducing an initial amount of syrup, enough to completely cover the tubes of the heat exchanger. Simultaneously vacuum pressure is established and the heating process is started by activating the corresponding loops and establishing the adequate set points. Then, the concentration phase follows with the objective of reaching the required degree of concentration, as reported by an online sensor. The concentration to achieve is assumed to be in correspondence with some appropriate supersaturation in the metastable region ( Figure 4D). As soon as the preconfigured concentration is reached, the sugar crystal seed is introduced. The seed consists of the population of sugar crystals that are going to be grown along the cycle. The next growing of the grain stage is the lengthier and most important of the cycle. Its task is to keep the supersaturation in the metastable region. Supersaturation is controlled in this stage by trying to strike the right balance between the rate of water evaporation and the feeding of liquor with the intention to make up for the sucrose that has migrated to the faces of the growing crystals. Because supersaturation is not directly measured, the required control is achieved by stipulating the concentration that the massecuite inside the chamber has to attain at every moment in the evolution of the cycle. This is done by establishing a concentration (brix) vs level curve that should be adjusted to each unit taking into account the expected purity of the feeding liquor. The growing phase finishes when the allotted space for the mass in the chamber is fully occupied. Then, a further concentration step follows. Finally, after breaking the vacuum, the mass is discharged into the common buffer tanks in its way to the centrifuges.
In Figure 5, the profiles of key crystallizer variables and parameters are shown, such as the level attained by the mass (Figure 5A), and the massecuite and mother liquor brix and purity in Figures 5C and 5D, respectively. Note the radical reduction of the effective heat transfer coefficient (U in W/m 2 K) of the heat exchanger ( Figure 5B) caused by a combination of factors to be discussed later on.  The length of the crystallization cycle is influenced by the concentration and purity of the incoming syrup as shown in Figure 6: Lower concentrations or purities lead to larger batch times.

Simulation model
The simulated case study offers a simplified version of the mathematical models, which have been taken from a full-fledged dynamical simulation developed for the operator's training. 38,40,41 They are, however, faithful enough to the actual plant behavior, especially in what concerns to the interface between evaporators and crystallizers. The description of the main circulating streams is a key element: There are the technical, impure syrups that are represented by three components: the solvent water and two solutes, ie, sucrose and impurities. The latter is treated as a pseudocomponent, described by a single parameter, namely, the purity of the stream, although in fact, it includes a wide variety of organic and inorganic substances. The massecuite stream adds the presence of sugar crystals, which are embedded in the syrup or mother liquor. The crystallization pan has been modeled as a homogeneous, perfectly mixed container. Dynamic balances to the whole mass (m) inside the chamber and for each of the components of the massecuite are carried out as represented in Equations (11)(12)(13)(14)(15). The mass flow rates (kg/s) entering and leaving the vessel are the feeding liquor (w feed ), the crystal seed (w seed ), and the massecuite discharged at the end of the cycle (w mc ).
The terms m suc , m imp , m cris , and m water stand for the masses of sucrose, impurities, sugar crystals, and water, respectively. The purity of incoming feed syrup and its brix are p ml and b ml , respectively. The first term makes reference to the relation of sucrose to all dissolved substances: sucrose plus impurities. The brix is the ratio of all solutes to water. On the other hand, p mc and b mc are the purity and concentration of the resulting massecuite, and in this case, the mass of crystallized sugar is considered, along with the dissolved sucrose, in the calculations. The crystal content term (cc) simply refers to the mass concentration of sugar crystals. Important ingredients of the balances are the mass flow rates representing crystallization (or dissolution) (w cris ) and the one accounting for the water evaporation (w evap ), one of the key factors involved in creating supersaturation. The term w cris , appearing in (12) and (14) and determined by using Equation (16), stands for the flow rate of sucrose leaving the solution to be incorporated to the surface (A c in m 2 ) of the growing population of crystals. The specific flux of crystallization (16) J cris (kg/s·m 2 ) is mainly dependent on supersaturation (s). In (17), supersaturation (s) is obtained as the ratio of the current mass concentration of sucrose (C suc ) to the sucrose concentration that marks its solubility in the technical solution (C suc_tech ). The latter is a function of the temperature (T ml in • C) and its purity (p ml ).
To model the dynamics of crystal growth, the known population balance equations (PBE) framework is typically used. The PBE is a partial differential equation depending on time and on a characteristic dimension (L in m) of the particle. However, in the case at hand, the exact shape of the CSD is not needed because only some aggregate values such as the surface of the population, used for example in (16), are required. Therefore, it is enough to follow the evolution of some moments of the size distribution by means of a set of ODEs as shown in (18). Only the first four moments have been considered. It is known that those moments have a direct physical interpretation: 0 (m −3 ) represents the number of particles, whereas 1 (m/m 3 ), 2 (m 2 /m 3 ), and 3 (m 3 /m 3 ) give us, respectively, the aggregate values of the length, of the surface and of the volume of the crystal population per unit volume of the solution. The terms Q seed and Q mc_mm (m 3 s −1 ) are the flow rates of the solution transporting the crystals during seeding and final discharge, respectively. The linear velocity of crystal growth is given by G (ms −1 ), whereas B (m −3 s −1 ) is the number new particles, of initial size L 0 (m), created each second per unit of solution volume. The term seed_k stands for the moments of the seed population, whereas V ml (m 3 ) is the volume of the solution. Dissolution effects are considered as a reduction of the size of the crystals, but the discrete event of their disappearance is not modeled. Other discrete phenomena such as agglomeration and breakage of crystals are likewise left out. The growing rate G is considered to be independent of size.
) + 1 · 10 7 · max(0, (s − 1)) 0.8 · cc 0.8 (20) There is an explicit relation between the flux of crystallization (J cris ) used in (16) and the linear growth velocity G, which is shown in (19) and is calculated taking into account the density of sugar crystals and two shape factors relating the characteristic size L with the crystal surface (f s ) and volume (f v ). The nucleation term (B) is calculated by means of (20), by adding the contributions of primary and secondary nucleation. Primary nucleation is the mechanism that works in the absence of other crystals, whereas secondary nucleation accounts for the presence of an already crystallized mass of sugar resulting on its current concentration (cc). The expression (20) has a form similar to the ones recommended in the work of Puel et al, 42 but the coefficients are chosen for the particular case at hand. The balance to the mass, which is evaporated from the boiling product and that shares the vessel volume with the latter could be described as in (21). The flow rate leaving the vessel is determined by a valve that is actioned by the control program of the crystallizer. The evaporated flow rate (W evap ) is pragmatically determined by means of (22), where P vap is the pressure of the mass of vapor and P eq is the steam equilibrium pressure for the mother liquor. The latter term takes into account the boiling point elevation at the given concentration and purity. The multiplying factor (K evap ) is arbitrary and should be selected just high enough to obtain the practical equality of both pressures. The vapor pressure at the top of the vessel is calculated using the known ideal gas law (23), knowing the mass of vapor, its temperature (T), expressed in K and assumed identical to the one of the boiling mass (T mc ). The volume V vap (m 3 ) can be determined from the algebraic expression (24) because the vapor competes with the massecuite occupied volume (V mc ) for the full capacity of the vessel (V). In (23), mm H20 (kg/mol) is the molar mass of water and R (J·mol −1 ·K −1 ) the ideal gas constant, whereas in (24), mc is the massecuite density (kg/m 3 ).
The energy balance to the product inside the chamber is given by means of (25), where h syp , h seed , and h mc are, respectively, the specific enthalpies (J/kg) of the streams of feed liquor, of the seed and of the resulting massecuite. The relationship between temperatures and enthalpies for these technical streams is tabulated for the given purity, brix, and crystal content. On the other hand, the term h evap represents the enthalpy of the evaporated water, which is a function of the existing pressure in the chamber (p vap ) and of massecuite temperature (T mc ). The determination of the heat transfer rate Q(W) (26) and the very related balance to the mass of steam in the primary of the heat exchanger (27) are both important elements of the energy balance (25) and of the modeling of the steam consumption behavior of the crystallizers. The term Q results from the difference of temperature at both sides of the shell and tubes exchanger of total area A (m 2 ). The already mentioned temperature of the massecuite (T mc ) is obtained as a function of its enthalpy from (25), whereas the determination of T cal demands to carry out a balance to the heating steam mass (m cal ) in the primary of the exchanger as in (27). The flow rate of the heating steam (w st ) is determined, at one end, by the pressure delivered by the relevant evaporator effect, at the other by the pressure developed as a consequence of the accumulation of steam in the shell of the exchanger, and by the aperture of the controlling valve placed in the circuit connecting both. The heating steam would be considered saturated for simplicity although in practice it is clearly superheated because it is produced by evaporation in solutions that exhibit an appreciable elevation of the boiling point. The control of the valve, as it should be emphasized, is managed by a PID whose reference is decided by the crystallizer's control program. The flow rate of condensed water leaving the shell (w cond ) is calculated (28) as the ratio of Q and the specific latent heat of condensation/evaporation h cond (J/kg), which is an exclusive function of the pressure p cal if the saturated steam assumption is made. The value of p cal could be obtained from the gas state law, similar to (23), that could be considered valid also for this case because the pressures involved are not too high. For the case at hand, the volume of the exchanger shell would be a fixed parameter, and the mass would be taken from (27) to then determine P cal . Note that in this case, the temperature T cal would be readily obtained from saturated steam tables.
Special care has been given to reproducing the changing overall effective heat transfer coefficient U shown in Figure 5B. This behavior can be explained as the combined effect of multiple factors: The increase of the level of the mass implies a higher pressure in the vicinity of the exchanger and thus a reduced local temperature difference at both sides of its walls. It is also the case that, as the viscosity increases and the evaporation induced bubbling is reduced, the mass circulation gets poorer causing a further decrease in the effective value of U. To model this effect in a very simple way, the expression (29) has been used. It uses a relation between U, the mass temperature, and its concentration adapted from the one proposed in the work of Semlali Aouragh Hassani et al 43 in relation to sugar evaporative continuous crystallizers. Here, however, an additional factor is added that decreases exponentially as the height of the mass (h) increases. The reduction starts acting over the midpoint of the total available height (h max ), approximately corresponding to where the upper rims of exchanger tubes are typically located, which is considered to be the level of optimum mass circulation.
The complexity of the model, even in a simplified version, is considerable. The task is, however, manageable using object-oriented modeling and simulation techniques that allows for the interconnection of parameterized instances of models by means of compatible ports in what is known as object-oriented (OO) or physical modeling. 44 This approach considers the hierarchical construction of models, reusing previously developed compatible objects such as valves, tubes, PID controllers, and tanks among other ancillary equipment. The circuitry conformed by product and steam conducting tubes, valves, and other connecting elements were not explicitly described above but are automatically inserted as requested by the described plant topology. Similarly, the parallel battery of crystallizers is deployed as three instances of the same basic model, appropriately connected.
The OO tool is able to deal directly with mathematical models expressed in the form of differential algebraic equations (DAE). In preparation to the numerical integration, the automatic determination of the computational causality of the model is performed. The appropriate ordering of equations is found, along with the identification of all the algebraic boxes that might exist and that would require an iterative solution at each integration step.
The evaporator station follows the same principles of physical modeling composition. The key element, the individual evaporator effect, has been developed following assumptions similar to those already described concerning the evaporation in crystallizers. Evaporator effects, however, work in a continuous fashion, regulating the level of mass inside the chamber by controlling the amount of juice leaving the unit. The required mass and energy balances are performed and the overall heat transfer coefficient is kept constant. Each effect steam output is connected to the primary of the following effect exchanger. The last effect (III) evaporated steam is delivered to a barometric condenser via a control valve that is driven by the slave PID of the cascade combination used to regulate the concentration of the syrup delivered at the output. The effects I and II provide the steam not only to the following evaporation effect but also to other consumers in the factory, but mainly this is a key element of the proposed plantwide problem to the crystallizers.
An indispensable ingredient of any process industry model is the use of physical and chemical properties for the different streams involved. The dependence of supersaturation on purity and temperature for technical sucrose solutions is important. Properties such as densities, specific enthalpies, boiling point elevation, among many others, as applicable to syrups, massecuite, sugar crystals, and steam, are also required. The compilation offered in the work of Bubnik et al 45 has been used as the main source for this information.

OPTIMAL PLANTWIDE PROBLEM OF THE SUGAR CASE STUDY FOR ENERGY EFFICIENCY
In the proposed sugar case study, the fed-batch parallel crystallizers are contending for the resources represented by the syrup to process, kept provisionally in the liquor buffer tank, and for the capacity of the corresponding effect for providing steam at the required pressure. The complexity of each individual crystallizer justifies the use of an optimizing controller for its optimal management. The design and implementation of this kind of controller is not a trivial task. There are contributions exclusively dedicated to the control of batch crystallizers in general and of sugar ones in particular. 46 Here, for the sake of simplicity, we assume that the existing recipe-based control, based on the brix versus level curve, is good enough. However, to give some control over the pace of crystallization, we introduce a factor P Fi centered at one that serves the purpose of multiplying the recipe-proposed set point for the pressure of the heating steam that comes from the evaporator effect to each one of the crystallizers. Note in Figure 7A, the effect of an increasing P F in the form of a displacement, parallel to itself, of the pressure set point at the primary of the heat exchanger, and as a result of it, the acceleration or slowing down of the pace of evaporation, and the resulting modification of the length of the cycle. Figure 7 also shows the impact of the evaporation pace in the supersaturation ( Figure 7B) and on the profiles of the flow rates of the streams that interact with the common shared resources: the feed liquor ( Figure 7C) and heating steam flow rate ( Figure 7D). Note the very pronounced time-varying characteristics of mentioned profiles. Discounting the initial syrup charge, both flow rates show a similar decaying pattern decided mainly by the already described reduction of U and reflecting the discussed supersaturation control scheme.
Working at a slower pace, with lower steam pressures and flow rates along the cycle, is economically beneficial in the sense that it reduces steam consumption. The multiple ways in which the previous assertion is valid will be discussed later on.
Moreover, a crystallization performed at lower evaporation rates implies working with reduced supersaturations. This results in a better CSD of the sugar product with favorable repercussions in their perceived market quality. Summing up, the task at hand is to schedule the cycles of the semibatch crystallizers units and the pace of the evaporation (decision variables T loadi and P Fi ) to meet two possible scenarios as follows.
1. Partial_plant case: improving steam efficiency by going at the slowest evaporation rate that is still compatible with the processing of the arriving syrup, which is decided by the upstream part of the factory. 2. Full_plant case: Add to the previously mentioned decision variables, the flow rate of the incoming juice, to state the optimal operation as the one that increases the plant throughput until a limiting restriction related to some unit capacity (a bottleneck) is met. This can also be done in the most efficient manner, using only the strictly required amount of steam energy.
Both scenarios can be stated in the proposed two-layer coordinated framework developed in Section 2 ( Figure 8). In what follows, however, the partial plant variant of the problem has been chosen. Note that it corresponds better to the real position of the subplant in the flow sheet of the whole sugar plant. The brix set point (bxSp) and the pressure (p vap_in ) at the input of the evaporation station, respectively (Figure 3), could have been taken as additional degrees of freedom. It has been decided, however, to keep them both fixed. In particular, the brix set point has been put to a high safe value as

Sugar case study two-level coordinated plantwide optimization
The optimization problem of the coordinator, as applied to the particular case at hand, is represented by Equations (30)(31)(32)(33)(34)(35)(36)(37)(38)(39). The performance function is considered to be the aggregate value of the consumption of steam by all the crystallizers in the whole exercise. Each crystallizer is considered to perform an identical number of cycles. There are nc crystallizers units, indexed by i, and nb cycles for each unit with index j. Note that it would have been perfectly defensible to seek for the minimization of another variable, such as, for example, the steam flow rate entering the evaporation station (W vap_in ). min.
The decision variables are the initiation instants for each one of the batches into the future for each crystallizer (T loadij ) and the coordinating variables between the upper and the subordinate controllers. The latter has been taken here, for convenience of implementation, as the maximum allowed batch time conceded to each subordinate unit (T b_mxij ). Note that there is a direct monotonic relationship between weighting the heating steam flow rate with a higher relative cost, as represented by the corresponding normalized parameter i in (8), which promotes working at a lower average value of that flow rate along the cycle, and the resulting slower pace of evaporation, which corresponds to a larger allowed maximum batch time (T b_mxij ). In (31)(32)(33)(34)(35)(36)(37), a very simplified surrogate model of the common inventories and of the relevant behavior of the consumption of those common resources from the crystallizers acts as a dynamic equality constraint for the problem. The dynamic model of the normalized level (L buf ) in the feed syrup buffer (31) is determined from a balance of the flow rates of the incoming concentrated syrup (w feed_conc ) and the feed syrup consumption of the installed crystallizers (w feed ). In (31), A buf represents the area of the base of the tank and L buf_max is maximum height. The syrup density feed has been taken, for simplicity, as the average value, which is approximately valid for the range of concentrations and purities expected. The inventory of the other shared resource, heating steam, is characterized by the pressure of the vapor space of the concerned evaporator effect described as a function of the drawn steam flow rate. A detailed account of this behavior would require the modeling of the complete effect with mass and energy balances in the syrup and vapor spaces. This more complete model would be the way to go if other decision variables, directly related to the evaporation, eg, p vap_in or bxSP, or the full-plant case, were to be considered. Here, however, because the interest lies only in the representation of the limited power of the steam source, the very simple static expression given by (32) has been used. It states that the instantaneous pressure provided by the heating steam source at any given time (p st ) falls exponentially as a function of the total combined instantaneous flow rate demanded by the crystallizers (w st ). It goes from a value established by the parameter p st_max , in such a way that, when w st reaches the parameter w st_max , the pressure has dropped all the way down to zero. The surrogate model not only needs to have access to the current estimation of the profiles of the of syrup (w feed ) and steam (w st ) consumptions, but it should also be able to generate dependable predictions of the variation of those profiles as a function of decision variables of the upper-level optimization: T loadij and T b_mxij . The option of using duplicates of the low-level crystallizer models to accurately reproduce the aggregate behavior of the crystallization section is out of the question. That solution would defeat the intention of having a parsimonious, noncentralized, loosely coupled two-level solution. Then, the alternative here proposed, already presented in a more general setup in Section 2.1, is based on obtaining, from each one of the subordinate controllers, the most recent prediction of the profiles of syrup (W feed_aij ) and steam (W st_aij ) consumption and their sensitivities to their respective coordinating variable for each unit and cycle ( W feed / T b_mx and W st / T b_mx ). The mentioned profiles are obtained from the previous value of the coordinating variable T b_ma . With this information, the linear approximations represented by (33)(34) should be able to predict how each of the profiles (W feedij , W stij ) would be modified for small variations in the corresponding coordinating variable T b_mx . It should be noted that both Equations (33)(34) are dependent on the local variables τ ij , which express the elapsed time within each cycle. They are defined for local times taking values from zero to the maximum cycle time of each crystallizer (T maxi ). The inequality restrictions imposed by Equation (35) establish a trust region for the validity of the linear approximations, allowing excursions of the coordinating decision variables in a range of 5% centered at the previous values (T bm_ aij ). The two fixed parameters T maxi and T mini serve the purpose of providing absolute maximum and minimum limit values for the range of possible batch cycle durations for each of the units.
The actual predicted flow rates of syrup (w feed ) and steam (w st ), as needed in (31)(32), should be defined in the real global time t for the complete horizon. These predictions are built by means of Equations (36) and (37), respectively. Both profiles are assembled by aggregating displaced versions of the proposed local profiles for all batches in all crystallizers. The displacements are provided by the decision variables T loadij currently established in the search for the optimal scheduling of the fed-batch subordinate units. Note that, for Equation (36)(37) to work, the value of the local profiles (W feedij , W stij ) in Equation (33)(34) should be considered as zero outside the range [0, T maxi ].
Finally, the inequality constraints (38) and (39) are required to make that all the decisions taken are compatible with keeping the level in the syrup buffer tank in the range [L min , L max ] and the pressure of steam in the evaporator effect above the value P st_min . The election of p st_min should be above the higher value of the set point of the heating steam pressure PID loop proposed by any of the crystallizer's program, taking into consideration the higher possible multiplying acceleration factor (P Fi ). If constraints (38)(39) are met, the subordinate problems can be solved independently.
The coordinating variable is put as a hard constraint on the allowed cycle duration by means of (43).

Sugar case study two-level coordinated plantwide optimization
The discussed plantwide optimizing control is going to be implemented by EMPC controllers at the two levels. The horizon of the coordinator EMPC is established, not as a fixed time, but as a given identical number of batches to be processed by each crystallizer in the immediate future. On its turn, the subordinate EMPCs use a shrinking horizon scheme, typical of batch processes. The discussed loosely coupled optimization strategy allows the deployment of the controllers to different network-connected computers. Figure 8 shows graphically the implemented setup, highlighting the required interchange of information between the two levels.
The outline of the algorithms to be employed by the coordinator and by each of the subordinate controllers are given in Algorithm 1 and Algorithm 2, respectively. They implement a discretized version of their respective optimal control problem, using the known receding horizon strategy. Note that the sampling time used by the coordinator or by any instance of the subordinate controllers need not be same.
The required interchange of information between the two levels is always performed in a nonblocking, asynchronous fashion. The coordinator, for example, need to obtain from the subordinates, at every sample instant (indexed by k real ), the latest values of the profiles of common resources consumption flow rates and of their sensitivities to the current values of the coordinating variable. In addition, these need to be done for each crystallizer i and each batch j for the considered horizon. The request (lines 3-5 in Algorithm 1), in some cases, would take almost no time because the latest version of some of these profiles could be directly retrieved from a globally known repository implemented, for example, in a shared computer memory segment. Once the optimization problem is solved (line 7 Algorithm 1), the solution sequence for the scheduling and coordinating variables is also written in the common memory. The first elements of these sequences are the ones that are going to be applied for the current control as required by the receding horizon scheme. Similarly, each one of the subordinate controllers reads the most recent values for scheduling commands and coordinating variables from the common memory (line 2 in Algorithm 2) and, after the local optimization is carried out (line 12 in Algorithm 2), it writes the current profiles of the consumption flow rates and sensitivities to the same repository (line 13 Algorithm 2).
Each subordinate's EMPC (Algorithm 2) is also continuously executing in a loop that is triggered at a rate defined by its particular sampling time (indexed by k local ). The optimal control of the unit for the requested batch implies the solution of the problem given in Equations (40-43) (line 11 in Algorithm 2) after reading the most recent value for the coordinating variable (T b_mx ) from common memory (line 2 in Algorithm 2). Note that the controller could be already engaged, controlling the physical unit (local auxiliary variable in_batch is true). In that case, the within-batch optimal control is required, taking into account the current index k local . In this latter case, the first element of the optimal sequence should be applied for the control of the real underlying unit (line 14 of Algorithm 2).
However, subordinate controllers not only should be capable of dealing with the current batch cycle. They should also be able, when being in the idle state (in_batch in false), to start a new cycle when the global real time (t) is greater than the current relevant T load scheduling command. However, subordinates also need to answer requests concerning their predicted behavior when processing batches at arbitrary conditions into the future, to satisfy the coordinator's need of having dependable information inside the limits of its prediction horizon. Algorithm 2 can answer to all these different concurrent needs by keeping in a local repository the information concerning the state of the current execution (eg, in_batch, k local index variables) relevant to each particular instantiation.
In addition to the actions already explained, both types of controllers must, in general, obtain the values of key variables from field sensors to perhaps adjust their respective models or perform the model-based observation of the state of the controlled process (line 6 in Algorithm 1 and line 3 in Algorithm 2).
It should be emphasized that, with the discussed asynchronous scheme in place, all the agents involved can interchange information at their own pace, obtaining always the most recent version of the pertinent data. This allows the election of their respective sampling frequency with relative independence. It is true, however, that the design of the plantwide solution should not disregard the danger of using excessively outdated information. Therefore, this type of consideration should also be heed when choosing sampling frequencies or by designing any kind of synchronization mechanism on top of discussed algorithms. In this context, the substitution of fixed-time sampling with event-based techniques could be an interesting option.
A key element of the two-level strategy hinges on the determination and use of the sensitivities of the flow rates of shared resources with respect to the coordinating variables. They may possibly be obtained by having each local controller to extract the pertinent information by means of postoptimality analysis techniques. For the case at hand, however, the sensitivities were obtained by solving the local optimization twice, one time with requested value of the coordinating variable and the second time with a slightly modified one. The disadvantage of spending more computation time is compensated by the flexibility obtained and the possibility of solving the local problems in parallel.

DISCUSSION OF SIMULATION RESULTS
The proposed two-level coordinated plantwide solution, as applied to the sugar crystallization partial-plant case, has been put to the test in simulation. The results to be discussed correspond to an exercise that covers approximately one day. The four EMPC controllers for the coordinator and three parallel identical crystallizer units, as shown in Figure 8, have been implemented.
The coordinator's EMPC uses prediction and control horizon that are equal in length, consisting of two complete batches into the future for each one of the three crystallizers.
On the other hand, all the subordinate units' EMPCs use the shrinking horizon scheme: The prediction horizon always covers from the current sampling instant in the cycle until the end of the batch. The optimal control profile for the only decision variable that has been considered (P F ) is parametrized with a sequence of piecewise constant functions. The resulting control horizon has been chosen to consist of two stretches: The first has a length equal to the sampling time, whereas the second covers all the way to the end of the cycle.
As discussed previously, the sampling period for each of the controllers can be chosen independently to meet the dynamical requirements of the process under control. Therefore, for the case of the coordinator, due to the capacity of the available syrup buffer (200 m 3 ), the expected product flow rate, and the process capacity of the array of crystallizers (55 m 3 ), a sampling period of 30 minutes has been found to be adequate. On the contrary, the control of the batch crystallizers would, in general, demand a smaller sampling period. For this process, a sampling time of 2 minutes would meet the demands dictated by the dynamics of the batch units. It is to be noted, however, that in this particular case, the subordinate EMPCs are not intended to control directly all the relevant details. The underlying low-level recipe-based control program assumes that responsibility, and the EMPC is only in charge of steering the unit toward a working regime of low-steam flow rate consumption, resulting in the longest possible cycle duration as demanded by the upper-level coordinator. Therefore, here, nothing is really gained by using a higher sampling rate because the changes of the subordinates control variables (P Fi ) end up as being naturally synchronized at the frequency of the requests made by the coordinator. Figure 9A shows the mass flow rate of juice entering the evaporation station, whereas in Figure 9B, its purity, also a known perturbation, makes the excursions shown in Figure 9B. It should be recognized that both inputs are subjected to simulated changes that are much more abrupt that those to be expected in the real situation. Other significant input variables have been kept fixed, eg, the juice brix at 22%, its temperature at 123 • C, and the set point of the brix controller at the output of evaporation (BxSp) at upper limit 76% of the typical range of [67%-76%]. The pressure of live steam pressure at the input of evaporation (P vap_in ) with typical values in the range [0.19-0. 22 MPa] has been fixed at 0.21 MPa. The maximum (T mx ) and minimum (T min ) cycle times should be chosen individually, taking into account characteristics of the specific unit for this case they have fixed at 12 000 and 7000 seconds, respectively, for all the batch units.
In Figure 10A, the decision variables for the scheduling, expressed as the predicted time for giving the starting commands to each of the subordinate crystallizers, and measured since the beginning of the exercise, are recorded. Figure 10B, on its part, gives the evolution of P F factors, the decision variables of each of the subordinate EMPCs. Note that all the P F factors are, for the most part, at the lower limit. This is a result that is consistent with the objective of working at a slower pace with the minimum consumption of steam. In the present optimization exercise, the P Fi decision variables are allowed to change, from one sampling time to the next, to any value in the allowed range [0.7-1.3]. In some applications, considerations regarding the convenience of a smoother operation of the crystallizers would advise to limit the acceptable changes to a more reduced neighborhood of the current value. This would obviously require the addition of the necessary restrictions to the optimization problem.
In Figure 11A, the evolution of the simulated value of each crystallizer supersaturation shows that it complies, during the periods of seeding and growing of the grain, with the metastable zone limits, approximately in the range [1-1.35] that  assures a correct CSD and, consequently, a good product quality. In Figure 11B the evolution of the level of the mass inside each pan gives a more intuitive assessment of the proposed scheduling.
The restrictions of having the level of the buffer tank ( Figure 12A) always inside the range [10%-80%] and the pressure of the heating steam source above 0.16 MPa ( Figure 12B) implies that it is legitimate to treat the predictions of the parallel subordinated EMPC as noninteracting. In the case of the steam pressure restriction, and to be on the safe side, in the simple model for the steam source used in the coordinator (32), the parameters p st_max and w st_max have been taken as 0.21 MPa and 20 kg/s, respectively, with the intention of underestimating the power delivering capacity of the real evaporator effect. Finally, in Figure 13, the coordination variable between the upper and lower layers is shown for the relation between the coordinator and the first subordinate crystallizer. It gives the evolution of the restriction on the subordinate EMPC-allowed cycle length. Note that, for almost all cases, the upper limit of the allowed trust region range for the linear approximation validity of the surrogate model is an active constraint of the optimization problem. The dynamic optimization has been approached by means of the direct sequential approach or control vector parametrization. 47 Only the control variables have been discretized, whereas the differential and algebraic equations are integrated numerically along several iterations in search of the locally optimal solution. The models representing the real plant, and those used by the controllers and the logic of the described two-layer strategy have been implemented using the already mentioned EcosimPro/Proosis tool. 48 The resulting nonlinear mathematical programming problems have been solved using the sequential quadratic programming-based numerical algorithm SNOPT. 49 The very core of the described scheme demands the use of surrogate models in the coordinator, which necessarily differs from those representing the real plant. In the case of the subordinate EMPCs, however, no mismatch between models and plants has been considered.
The average solution time for the coordinator problem is of 312 seconds each sampling time. The subordinate EMPC consume much less time, around 10 seconds, for the simple, one-decision variable here presented. The results were obtained with an Intel CORE i7-powered PC, with a clock frequency of 1.99 GHz and 16 GB of RAM. In what concerns the scalability of the plantwide controller, the possibility of distributing the computational burden should be taken into account, which is designed in the very loosely coupled nature of the proposed strategy.

Energy saving considerations
In most sugar factories, the scheduling of crystallizers is performed with the direct intervention of control room operators. They need to keep track of the flow rates of the relevant streams, the level of the buffer tanks, and the current stage in the evolution of each crystallizers' cycle, among other factors. With this information, they are able to issue the necessary commands to each unit, including the modification of the relevant parameters of the control programs and the batch cycle starting instructions. In this context, to keep the problem manageable, operators tend to conduct the crystallizers at the higher pace that is still compatible with the availability of syrup and steam. Therefore, to discuss the energy-saving implications of the two-layer EMPC-based solution, it seems fair to compare the results obtained under the proposed strategy with those gathered from simulations carried out for the same conditions but with different fixed evaporation rates. In Table 1, the data obtained for the optimal controller sequence of P F values, shown in Figure 10B, are registered (as P F opt.) alongside the results of three simulations performed at nominal conditions (P F = 1), with maximum evaporation rhythm (P F = 1.3) and with the minimum possible one (P F = 0.7). It is to be observed that the EMPC control outperforms both the nominal and fast-paced regimes, in the sense that it requires a lower aggregate total consumption of steam (in T) in conformity with the cost indexes of all the solved optimization problems. It is also better in what concerns the total heating energy demanded (in MJ), and in the average ratio of that energy to the total amount of syrup processed in the crystallizers (MJ/T).
On the other hand, the results exhibited by the EMPC plantwide solution compares very positively with the lower possible steam consumption attainable with a P F control variable in the inferior limit of the allowed range (P F = 0.7). The average steam energy per unit of syrup mass processed by all crystallizers, for example, is only 0.40% higher than the one of the most favorable case, but almost 4.1% better than those corresponding to the higher rates. It should be noted that all the simulations taken for comparison have been obtained using the scheduling proposed by the implemented controller, but considering a limitless supply capacity of syrup and steam. Table 1 also shows that a lower energy demand from the crystallizer is translated as a correspondingly reduced steam intake at the evaporators input. It should be said that to increase the realism, the model assumes the presence of two streams of vapors, produced by effects I and II of the evaporation, that serves the purpose of simulating the demand from other nonmodeled units of the factory. Therefore, to interpret the data shown correctly, it should be considered that, in addition to the steam sent to the crystallizers, there are these two other steam flow rates that have been made proportional (9% and 7% respectively) to the rate of juice to concentrate entering the department. However, the benefits of conducting the crystallizers at lower heating steam pressures (lower P F factors) are not limited to the results just discussed. To take full advantage of the energy reutilization scheme of sugar multieffect evaporation, it is recommended to take the heating steam from the effect with the lower pressure, which is still compatible with the heating needs of the consumer. In the case study considered, this would imply that it is much better to use as steam source effect II, as compared to the original effect I (Figure 3). Of course, this change is only possible if the temperature and pressure requirements of the primary of the crystallizer's heat exchangers are low enough as to be met by the new lower-temperature steam source. This condition is satisfied in the benchmark and by many modern sugar installations. Therefore, an important advantage of working at low pressures and evaporations rates would be that it opens the door to improve the plant's economy by making the mentioned steam source re-assignations. Table 2 shows the impact of changing the heating steam source. Both simulations were performed with the scheduling and P F factors provided by the EMPC controller. The evaporation effectiveness rate (K) makes reference to the ratio of the aggregate energy of the vapors supplied by all the evaporator effects in relation to the amount of live steam energy required by the first effect. It gives a measure of how many times the original energy from the boilers is reutilized, so a larger number, as the one obtained in the case at hand (as shown by corresponding row in Table 2), is always better. The improvement of using effect II is revealed also in the reduction of the total mass of live steam consumed by the evaporators for the complete simulation exercise. Likewise, the average specific value of live steam consumption per mass of syrup processed (in MJ/T) is also notably improved by the re-assignation.
An in-depth treatment of the energy management of a complete sugar factory is beyond the scope of this article. However, in general, it can be said that the efforts for improving energy efficiency should strive for a better reutilization in the evaporators (K parameters as high as 3.5 are possible with five effect evaporators), but they also involve many other issues, including a carefully designed process integration, incorporating the recycling of the different condensate streams. 50 In this context, it is considered that, due to the low temperature involved, the steam sent from the last effect to the condenser is completely wasted. It is for this reason that a low enthalpy associated with this stream would be a good marker of the factory's efficiency. The simple reduction of steam consumption from the crystallizers does not affect this indicator in a significant way, but the steam source re-assignation, made possible by that same reduction, does. This is made explicit in Table 2, showing the considerable decrease of effect III vapor mass flow rate to the condenser and of its specific energy per mass of concentrated syrup.
On the other hand, the use of lower-pressure effects is also favorable from an exergy balance perspective. It liberates capacity of the higher-pressure effects, an energy that can then be used on the heating of higher temperature streams.
Finally, the existence of many important challenges concerning the efficient use of energy in the sugar industry should be recognized. For example, in most modern factories, the steam from the boilers is not used exclusively for heating purposes. It is also employed to power back-pressure turbines in what constitutes an electricity cogeneration system. The electrical energy obtained with this arrangement is frequently able to supply the needs of the factory, and in cases of surplus, it can be sold for a profit, to the external grid. 51 The optimal harmonization of this cogeneration capacity with the heating needs of the crystallizers and of the rest of the factory is a promising open problem and represents a natural continuation of the work here described.

CONCLUSIONS
A novel plantwide hierarchically coordinated economic optimization for process plants with continuous and semibatch processing units sharing common resources has been proposed. The methodology has been described in the context of energy efficiency and applied to a realistic beet sugar plant simulation. The nature of the proposed solution would allow, with certain ease, the dynamic addition or removal of the subordinate units, which are loosely coupled, in a plug and play style, which is in tune with what is currently advocated by the proponents of the smart factory concept.