Second-order cone optimization of the gradostat
aa r X i v : . [ m a t h . O C ] F e b Second-order cone optimization of the gradostat ⋆ Josh A. Taylor
The Edward S. Rogers Sr. Department of Electrical and Computer Engineering, University ofToronto, Toronto, Canada
Alain Rapaport
MISTEA, Univ. Montpellier, INRAE, Institut Agro, Montpellier, France
Abstract
We maximize the production of biogas in a gradostat at steady state. The physicaldecision variables are the water, substrate, and biomass entering each tank and the flowsthrough the interconnecting pipes. Our main technical focus is the nonconvex constraintdescribing microbial growth. We formulate a relaxation and prove that it is exact whenthe gradostat is outflow connected, its system matrix is irreducible, and the growth ratesatisfies a simple condition. The relaxation has second-order cone representations forthe Monod and Contois growth rates. We extend the steady state models to the case ofmultiple time periods by replacing the derivatives with numerical approximations insteadof setting them to zero. The resulting optimizations are second-order cone programs,which can be solved at large scales using standard industrial software.
Keywords:
Gradostat; second-order cone programming; convex relaxation; wastewatertreatment; biogas.
1. Introduction
The gradostat is a nonlinear dynamical system in which multiple chemostats areinterconnected by mass flows and diffusion. In each chemostat, microbial growth convertsa substrate to biomass. This also produces biogas, a useful energy source. Our primarymotivation for this setup is the design and operation of a network of wastewater treatmentplants [1].We seek to maximize the production of biogas in gradostats with two standard growthrates: Monod [2] and Contois [3]. There are several physical decision variables, includingthe inflows of water, substrate, and biomass at each tank, and the installation of pipesbetween tanks. We start with a steady state model obtained by setting the derivatives ⋆ Funding is acknowledged from the Natural Sciences and Engineering Research Council of Canadaand the French LabEx NUMEV (Project ANR-10 LABX-20), incorporated into the I-Site MUSE, whichpartially funded the sabbatical of J. Taylor at MISTEA lab.
Email addresses: [email protected] (Josh A. Taylor), [email protected] (AlainRapaport) n the gradostat to zero. The resulting algebraic equations specify a feasible set, whichis nonconvex due to the nonlinear growth in each tank and, in some setups, discrete andbilinear mass flows between the tanks.Nonconvex optimization can be difficult even at small scales. Second-order coneprogramming (SOCP) is a tractable, convex optimization class that is often used toapproximate nonconvex problems [4]. In this paper, we construct SOC and mixed-integer (MI)SOC relaxations of the gradostat. Our main original contributions, listedbelow, center on the convexification of the growth rate constraint. • In Section 3, we formulate a simple relaxation of the gradostat. The relaxationis obtained by allowing the conversion of substrate to biomass to be less than orequal to the growth kinetics. In Theorem 1, we prove that when the gradostat isoutflow connected, its system matrix is irreducible, and a simple condition on thegrowth rate is satisfied, this relaxation is exact, which is to say that the inequalityis satisfied with equality. • In Section 4, we identify original SOC representations of the relaxed growth con-straints. Specifically, the Contois growth constraint is exactly representable as anSOC constraint. The Monod growth constraint is SOC under a constant biomassapproximation. We use the convex envelopes of [5] to obtain an SOC outer approx-imation of the Monod growth constraint in the general case. • In Section 6, we give two extensions. In Section 6.1, we give simple linear under-estimators of the growth constraints for when the relaxations are not exact. InSection 6.2, instead of setting the derivatives in the gradostat to zero, we replacethem with linear numerical approximations. This leads to an optimization withmultiple time periods, which can accommodate transient conditions.The end result is a family of SOCPs and MISOCPs for optimizing the gradostat insteady state. Today, SOCPs with tens of thousands of variables and constraints can besolved in seconds on a typical personal computer. MISOCPs are also reasonably tractablebecause, like mixed-integer linear programs, there are powerful mathematical tools forspeeding up their solution [6, 7, 8], and they are handled by industrial solvers such asGurobi [9]. This enables us to solve each MISOCP to optimality at moderate scales,typically up to one hundred binary variables.We now review some relevant literature. We refer the reader to [10, 11] for compre-hensive coverage of the chemostat. The gradostat was originally formulated in [12] as asingle series of tanks and later generalized to a network of interconnected chemostats.There is an ongoing literature stream focusing on its nonlinear analysis [13] and con-trol [14]. To date, there have been no applications of numerical convex optimization tothe gradostat in steady state.From a technical viewpoint, the nearest topic to ours is the optimization of chemicalprocess networks [15]. Common features include bilinear mass flows, which we similarlylinearize using disjunctive programming [16, 17], and quotients of variables, which weapproximate with convex envelopes [5, 18, 19] in Section 4.2.2. The main feature of thegradostat that is not present in chemical process networks is the microbial growth in thetanks, which is the core focus of this paper.The design of chemical reactors has been studied for decades, see, e.g., [20]. Ahandful of papers within this stream have focused on optimal design, some using the2hemostat and gradostat. An early text is [21], which uses dynamic programming, butnot the chemostat. Several later studies model growth with the Monod equation andderive analytical expressions for parameters such as concentration, residence time, andtank volume [22, 23, 24, 25]. References [26, 27, 28] build on this approach, derivingnumerical and qualitative conditions for optimizing reactors modeled as ‘steady-stateequivalent biological systems.’More recent studies have explicitly optimized the design and operation of intercon-nected tanks. The volumes of series bioreactors with Monod and Contois growth areoptimized in [29, 30] for a given output substrate concentration at steady state. The op-eration of two chemostats in series is optimized in [31], and the volumes and diffusion rateof a chemostat with a side compartment are optimized in [32]. Of particular relevanceare [33, 34], which model interconnected wastewater treatments plants as a gradostat.The latter formulates inflow management as a model predictive control problem, whichit solves using particle swarm optimization. While closely related to our perspective,these papers focus on different problem statements with specific network structures, anddo not employ convex relaxations or SOCP.The rest of the paper is organized as follows. Section 2 covers the relevant backgroundand states the nonconvex optimization problem. Section 3 gives a simple relaxationand analyzes when it is exact. SOC representations of the Monod and Contois growthconstraints are given in Section 4, and the bilinear mass balance constraints are linearizedin Section 5. Section 6 gives linear underestimators of the growth constraint and extendsthe models to the case of multiple time periods. In Section 7, we summarize the modelsand implement those that are SOCPs and MISOCPs in numerical examples.
2. Setup
We consider a gradostat with n interconnected tanks, and denote the set of tanks N .Tank i has water inflow Q in i and outflow Q out i . The concentrations of the substrate andbiomass inflows are S in i and X in i . The tank is assumed to be perfectly mixed and hassubstrate and biomass concentrations S i and X i . We suppress the subscript to denotethe vectors of these quantities in R n . V ∈ R n × n is a diagonal matrix in which V ii is thevolume of tank i .The substrate in tank i is converted to biomass at the rate r ( S i , X i ) /y = µ ( S i , X i ) X i /y ,where the growth rate, µ ( s, x ), is positive for s > x >
0. We refer to µ ( s, x ) x asthe kinetics. The conversion occurs with yield y >
0. We will focus on the following twogrowth rates. • Monod [2]: µ M ( s ) = µ max sK + s . • Contois [3]: µ C ( s, x ) = µ max sKx + s . We let all tanks have the same µ max , K , and y , and note that our results straightforwardlyextend to the case where they are not identical due to, for instance, temperature or pH3ariation. Some of our results apply to other growth rates as well. We refer the readerto Appendix 1 of [35] for a comprehensive list.The Monod growth rate, µ M ( s ), is concave. The corresponding kinetics, µ M ( s ) x ,are quasiconcave, but not concave. The Contois growth rate, µ C ( s, x ), is not concave.Observe that the corresponding kinetics can be written µ C ( s, x ) x = µ M ( s/x ) x . This isthe perspective of the Monod growth rate. The kinetics corresponding to the Contoisgrowth rate are therefore concave because the perspective of a concave function is alwaysconcave (see, e.g., Section 3.2.6 of [36]).We let Q ij denote the flow from tank i to tank j . Flow conservation implies that Q in i + X j ∈N Q ji = Q out i + X j ∈N Q ij (1)for i ∈ N . Let ˜ d ij denote the diffusion between tanks i and j , where ˜ d ij = ˜ d ji . Notethat this could represent the sum of the diffusions in multiple pipes. For example, if d ij is the diffusion in a pipe with flow from i to j , and d ji for another with flow from j to i , then ˜ d ij = d ij + d ji . We encounter this scenario in Section 5, in which one has thedecision to build a pipe in either direction.Let C = diag (cid:2) Q in i (cid:3) , G = diag [ Q out i ], and M ij = (cid:26) Q ji , i = j − Q out i − P k ∈N Q ik , i = j = (cid:26) Q ji , i = j − Q in i − P k ∈N Q ki , i = jL ij = (cid:26) ˜ d ij , i = j − P k ∈N ˜ d ik , i = j . Let ∈ R n be the vector of all ones. In matrix form (1) is given by ( M + C ) = 0.Observe that because L = 0, we also have ( M + L + C ) = 0. Similarly, ⊤ ( M + G ) = 0and ⊤ ( M + L + G ) = 0.The gradostat is a type of compartmental system, and M is a compartmental matrix.A compartmental system is outflow connected if there is a directed path from every tankto a tank with outflow, i.e., a tank i with Q out i >
0. A key property of compartmentalmatrices is that the matrix M is invertible if and only if the system is outflow connected. M is irreducible if it cannot be made block lower triangular by reordering its indices. Inthe gradostat, this means that there is a directed path from each tank to every othertank. Note that if the gradostat is outflow connected and M is irreducible, then M + L is respectively invertible and irreducible. We refer the reader to [37] for a thoroughdiscussion of compartmental systems.We henceforth assume that all gradostats are outflow connected. The matrices M and M + L are invertible, and therefore − ( M + L ) − C = (2a) − (cid:0) M ⊤ + L (cid:1) − G = . (2b)We also know that − M and − M − L are M -matrices [37]. This implies that − M − and − ( M + L ) − are nonnegative matrices. If M is irreducible, then − M − and − ( M + L ) − are positive matrices [38]. 4 .2. Equilibria of the gradostat The dynamics of the gradostat are V ˙ S = − y V R ( S, X ) + ( M + L ) S + CS in V ˙ X = V R ( S, X ) + ( M + L ) X + CX in , where R i ( S, X ) = r ( S i , X i ), i ∈ N . We obtain a steady state model by setting thederivatives to zero. We now briefly discuss when the solution to the steady state modelcorresponds to a unique, stable equilibrium of the gradostat.Let Z = X + yS . The dynamics of ( Z, X ) are V ˙ Z = ( M + L ) Z + C (cid:0) X in + yS in (cid:1) V ˙ X = V R (( Z − X ) /y, X ) + ( M + L ) X + CX in . Because the gradostat is outflow connected, M + L is negative definite and Z is globallyasymptotically stable with equilibrium ¯ Z = − ( M + L ) − C (cid:0) X in + yS in (cid:1) . Due to thecascade structure, the dynamics of X are asymptotically equivalent to those obtained byreplacing Z with ¯ Z . We assume that ¯ Z >
0; the following are two simple conditions thatguarantee this. • If C (cid:0) X in + yS in (cid:1) is not the zero vector and M is irreducible, then ( M + L ) − isstrictly negative [38] and ¯ Z > • Because M is a negative definite M -matrix, ( M + L ) − is negative definite andnonpositive [38]. If C (cid:0) X in + yS in (cid:1) >
0, i.e., there is inflow of substrate and/orbiomass at every tank, then ¯
Z > X in = 0, a ‘washout’ equilibriumwith X = 0 always exists. A positive equilibrium with S >
X > µ (cid:0) ¯ Z i /y, (cid:1) > − M ii − L ii for all i ∈ N . Intuitively, this means that the system canconvert substrate to biomass faster than it ejects biomass. In this case, the washoutequilibrium is repulsive. If r (cid:0) S i , ¯ Z i − yS i (cid:1) is increasing and strictly concave on (cid:2) , ¯ Z i /y (cid:3) ,then there is at most one positive equilibrium. It is easy to verify that this last conditionholds for Monod and Contois growth.The steady state approximation by definition limits the range of scenarios we canconsider. For example, a wastewater treatment system is clearly not in steady stateduring a storm surge. On the other hand, it may be a useful approximation when choosingwhere to install new pipes to improve efficiency under average operating conditions. Laterin Section 6.2 we drop the steady state assumption and instead numerically approximatethe derivatives. This enables us to optimize trajectories of the gradostat under time-varying conditions. The objective is to maximize the production of biogas. This amounts to maximizingthe conversion of substrate to biomass at a subset of output tanks,
M ⊆ N . Thecorresponding objective is max X i ∈M V ii r ( S i , X i ) . (3)5ll formulations in this paper are valid with the objective (3), but the theoretical resultsof Section 3 might not hold if M ⊂ N . If M = N , mass conservation implies that X i ∈N Q in i S in i = X i ∈N Q out i S i + 1 y V ii r ( S i , X i ) . If the left hand side is fixed and M = N , (3) is equivalent tomin X i ∈N Q out i S i . This corresponds to minimizing the substrate leaving the network.We hereon focus on maximizing the generic objective F ( T ), where T = R ( S, X ).This corresponds to biogas production if F ( T ) = P i ∈M V ii T i . It can also accommodateadditional features, e.g., concavity could reflect diminishing returns due to limited abilityto store biogas.
The full problem is given by P : max F ( T ) (4a)such that T i = r ( S i , X i ) , i ∈ N (4b)1 y V T = ( M + L ) S + CS in (4c) − V T = ( M + L ) X + CX in (4d)0 = ( M + C ) (4e) (cid:0) d, Q, Q in , S, S in , X, X in (cid:1) ∈ Ω . (4f)The variables are d , Q , Q in , S , S in , X , X in , and T . The growth constraint, (4b), equates T i to the kinetics in tank i . (4c), (4d), and (4e) balance the substrate, biomass, andwater in each tank. The set Ω in (4f) represents generic constraints on design and/oroperation. Several possibilities are listed in Use Case 1 below. Note that the tankvolumes are not variables, but could in principle be incorporated in a tractable mannerusing the techniques of Section 5. Use Case 1.
The following are operational constraints that could be represented by theset Ω . • X in = 0 . All biomass in the system either is converted from substrate or wasalready present before the system came to steady state. • For each i ∈ N , Q out i S i ≤ ˇ S i . The mass of the substrate released from each tankper unit of time cannot exceed some limit. • Q in ⊤ S in = ¯ S . The total mass of substrate that enters the network per unit of time, ¯ S , is allocated over the tanks. This could also be the case with biomass. ⊤ Q in ≤ ¯ Q . The total inflow (and, by conservation, ouflow) of water is limited.Observe that if the inflow is too small, the system will run inefficiently, but if it istoo large, washout will be a stable equilibrium. • For each i ∈ N , µ (cid:0) ¯ Z i /y, (cid:1) > − M ii − L ii . This ensures that the gradostat has apositive equilibrium and that the washout equilibrium, if it exists, is repulsive. Notethat here ¯ Z is an optimization variable, subject to the constraints M + L ) ¯ Z + C (cid:0) X in + yS in (cid:1) µ (cid:0) ¯ Z i /y, (cid:1) ≥ δ − M ii − L ii , where δ is a small positive constant that makes the inequality strict. If µ ( s, is concave in s , then this is a convex constraint. In the case of Contois growth, µ C (cid:0) ¯ Z i /y, (cid:1) is a constant. In the case of Monod growth, the latter constraint hasan SOC representation, which has the same form as that given in 4.2.1. The growth constraint in P , (4b), is nonconvex. Sections 3 and 4 construct SOCrelaxations of this constraint. We consider two cases for the rest of the problem. • When the variables d , Q , and Q in are constant, (4c) and (4d) are linear. Then theresulting relaxations of P are SOCPs. In this case, the physical decisions are theconcentrations of the substrate and biomass inflows, S in and X in . Note that S , X ,and T are also optimization variables, but are not under direct control of a systemoperator. • When d , Q , and Q in are discrete variables, (4c) and (4d) are bilinear. We usedisjunctive programming to linearize these constraints in Section 5. In this case,the relaxations of P are MISOCPs. The physical decisions are S in , X in , the flowsbetween tanks, Q , and the flows into the tanks, Q in . Q ij could represent, forexample, the decision to turn on a fixed speed pump in the pipe from tank i to j ,or the decision to build a pipe from i to j . We give more detail on the forms of d , Q , and Q in in this case in Section 5. Use Case 2.
Microbial activity in soil produces biogas, which, because it is not captured,contributes to the greenhouse effect. This is often modeled with Monod dynamics [40].It has also been shown to be dependent on spatial heterogeneity [41, 42], motivating theuse of compartmental modeling. While there are no true design variables, it is of interestto understand which spatial structures lead to the greatest release of biogas. We canestimate this by maximizing a gradostat’s biogas production over its water inputs, Q in ,and flows, Q and d . The biomass in soil evolves slowly [43] and is thus often treated asa constant parameter [41, 42]. We model this by dropping constraint (4d) from P andsetting X = X c . In Section 4.2.1, we show that in this case the kinetics correspondingto the Monod growth rate have an SOC representation.
3. Relaxation of the growth constraint
The growth constraint, (4b), is a nonlinear equality and hence nonconvex. We relax(4b) by replacing it with the inequality T i ≤ r ( S i , X i ) , i ∈ N . (5)7e refer to P with (5) instead of (4b) as P R . Unlike the equality (4b), (5) is convex forsome growth rates, and can sometimes be represented as an SOC constraint; this is thefocus of Section 4. In this section, we analyze when P R is exact, which is to say has thesame optimal solution as P .First, observe that if the optimal solution to P R satisfies (5) with equality for all i ∈ N , then P R is exact. Given a solution to P R , we can therefore determine if it isfeasible and optimal for P by simply checking if (5) binds. Definition 1.
Let
E ⊆ N be such that for each i ∈ E , Q ji + d ji = 0 for all j ∈ N , i.e.,it receives no flow from other tanks. A gradostat is fully fed if for each i ∈ E , Q in i > , S in i > , and X in i > . Lemma 1.
Suppose the gradostat is fully fed and outflow connected. Then
S > and X > . Proof.
Consider a path, L , from tank s to tank t . Assume that S in s > X in s >
0, and Q out t >
0, and that for each edge ij ∈ L , Q ij + d ij >
0. We proceed by induction on thepath.Observe that because flow must enter every tank and S in s > X in s >
0, we musthave S s > X s > s . Now suppose that jk ∈ L and that S j > X j >
0. Then S k > X k > k . Therefore, by induction, S i > X i > i ∈ L .Because we have assumed that the gradostat is outflow connected and fully fed, alltanks must lie on a path like L . Therefore, S >
X > (cid:3) In an outflow connected gradostat, there can only be zero substrate and biomass attanks which have no inflow of substrate or biomass and receive no flow or diffusion fromother tanks. Such tanks cannot exist in a fullly fed gradostat. Observe that anothersufficient condition for
S > L ) = n − Assumption 1. r ( s, x ) is concave, differentiable, and for all s ≥ and x ≥ , one has y ∂r ( s, x ) ∂s − ∂r ( s, x ) ∂x ≥ . We will discuss Assumption 1 for specific growth rates in Section 4.In Theorem 1 below, S , X , and T are the only variables in P R , the remaining variablesin P are regarded to be constant, and we drop the operational and design constraints in(4f). We will discuss how the theorem extends to the general case after the proof. Theorem 1. P R is exact if F ( T ) is concave and differentiable, the gradostat is outflowconnected, Assumption 1 holds, and either of the following is true: • M is irreducible and (cid:0) M ⊤ + L (cid:1) V − ∇F ( T ) ≤ and is not uniformly zero for all T ≥ ; or • the gradostat is fully fed and (cid:0) M ⊤ + L (cid:1) V − ∇F ( T ) < for all T ≥ . roof. P R is convex due to Assumption 1. If P R satisfies a constraint qualificationsuch as Slater’s condition, any optimal solution must satisfy the Karush-Kuhn-Tucker(KKT) conditions [36]. P R satisfies Slater’s condition if there is a feasible solution for which (5) is strict.Because the gradostat is outflow connected and either is fully fed or has an irreducible M , if there is any inflow of substrate and biomass, then we must have S >
X > r ( S i , X i ) > i ∈ N . We obtain afeasible solution for which (5) is strict by setting T = 0. Therefore, a Slater point exists,and any optimal solution of P R satisfies the KKT conditions.Let ρ ∈ R n + be the vector of dual multipliers of constraint (5), and let σ ∈ R n and ǫ ∈ R n be the respective multipliers of (4c) and (4d). Let U S and U X be diagonalmatrices with U Sii = ∂r ( S i , X i ) ∂S i , U Xii = ∂r ( S i , X i ) ∂X i for each i ∈ N . The KKT conditions for P R are given by(4 c ) , (4 d ) , (5) ρ = ∇F ( T ) + V (cid:18) y σ − ǫ (cid:19) (6a) U S ρ = (cid:0) M ⊤ + L (cid:1) σ (6b) U X ρ = (cid:0) M ⊤ + L (cid:1) ǫ (6c)0 = ρ i ( T i − r ( S i , X i )) , i ∈ N . (6d)The complementary slackness condition, (6d), implies that if ρ >
0, then constraint (5)binds for all i ∈ N and P R is exact.Let U = y U S − U X and W = (cid:0) M ⊤ + L (cid:1) V − − U . U is positive semidefinite due toAssumption 1. Because the gradostat is outflow connected, (cid:0) M ⊤ + L (cid:1) V − is negativedefinite, and therefore so is W . Arithmetic with (6a)-(6c) yields ρ = W − (cid:0) M ⊤ + L (cid:1) V − ∇F ( T ) , (7)Because − W is a positive definite M -matrix, W − is nonpositive. We have ρ > • If M is irreducible, W − is strictly negative [38]. Therefore, if (cid:0) M ⊤ + L (cid:1) V − ∇F ( T ) ≤ T ≥
0, then ρ > • If (cid:0) M ⊤ + L (cid:1) V − ∇F ( T ) < T ≥
0, then ρ > W − is negativedefinite and nonpositive.These are the two conditions we assumed in the theorem. Therefore, ρ >
0, constraint(5) is met with equality, and P R is exact. (cid:3) We now discuss Theorem 1 and its proof. We have regarded all variables except S , X , and T to be fixed. The theorem therefore indicates when a solution to P R , ( S, X, T ),solves equations (4b)-(4d). Theorem 1 also holds when d , Q , Q in , S in , and X in arevariables, so long as its assumptions hold at their optimal values. If constraint (4f)only affects d , Q , Q in , S in , and X in , it will not change the exactness. However, if (4f)constrains S or X , then Theorem 1 is not guaranteed to hold.9 orollary 1. Suppose that F ( T ) = X i ∈N V ii T i . P R is exact if the gradostat is outflow connected, Assumption 1 holds, and either of thefollowing is true: • M is irreducible; or • the gradostat is fully fed and Q out > . Proof.
The result is obtained by substituting ∇F ( T ) = V in (7) and using (2b) tosimplify the right hand side to ρ = W − Q out . (cid:3) There are likely further possible refinements and generalizations. For example, if F ( T ) is not differentiable, then the conditions of Theorem 1 need to hold for all of itssubgradients, or, alternatively, some subgradient at the optimal solution. On the otherhand, there are surely systems of interest for which the relaxation is not exact. From thispoint of view, Theorem 1 does not specify the entire set of gradostats for which P R isan exact relaxation. Rather, it is theoretical evidence that P R is exact for a meaningfulclass of gradostats and may be a high quality approximation for others.
4. SOC representations of P R In this section we construct original SOC representations and approximations of thegrowth constraint in P R , (5). With the Contois growth rate, constraint (5) takes the form T i ≤ µ max S i X i KX i + S i , i ∈ N . (8)This is a convex constraint because, as discussed in Section 2.1, the right hand side isconcave. Theorem 2. (8) is equivalent to the hyperbolic constraint ˜ S i + ˜ T i ≤ (cid:16) ˜ S i − ˜ T i (cid:17) (cid:16) µ max KX i + ˜ S i − ˜ T i (cid:17) (9a)0 ≤ ˜ S i − ˜ T i , (9b) where ˜ S i = µ max S i and ˜ T i = KT i , i ∈ N . This can be shown straightforwardly by simplifying (9a). Note that (9b) ensures thenonnegativity of both multiplicative terms on the right hand side of (9a), and is impliedby (8). We refer to P R with (9) in place of (5) as P RC .10yperbolic constraints are a special case of SOC constraints. In standard SOC form,(9a) is written (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ˜ S i ˜ T i µ max KX i (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ≤ µ max KX i + ˜ S i − ˜ T i . We now test Assumption 1 for P RC . We have1 y ∂r ( s, x ) ∂s − ∂r ( s, x ) ∂x = µ max (cid:0) Kx − ys (cid:1) y ( Kx + s ) . Clearly this is negative for some s ≥ x ≥
0, and hence does not perfectly satisfyAssumption 1. However, given that there is usually more biomass than substrate, y < K ≈
1, we expect it to hold around the optimal solution for most realisticsystems. Therefore, we expect Theorem 1 to hold most of the time for P RC . With the Monod growth rate, constraint (5) takes the form T i ≤ µ max S i X i K + S i , i ∈ N . (10)We refer to P R with (10) in place of (5) as P RM .Because the right side of (10) is not concave, Assumption 1 does not fully hold, andTheorem 1 does not directly apply. However, it is differentiable, which means that theKKT conditions hold at any local optimum. If the derivative condition in Assumption 1is satisfied along with the other conditions of Theorem 1 or Corollary 1, then P RM isexact. Evaluating the derivative condition in Assumption 1 yields1 y ∂r ( s, x ) ∂s − ∂r ( s, x ) ∂x = µ max (cid:0) Kx − ys ( K + s ) (cid:1) y ( K + s ) . As with the Contois growth rate, this is negative for some s ≥ x ≥
0, but is likelyto be nonnegative around the solution for a realistic system. Therefore, we expect P RM to be exact, albeit at a possibly local optimum.Unfortunately, without a more tractable representation of (10), P RM will be difficultto solve at larger scales. In Sections 4.2.1 and 4.2.2, we construct two different SOCapproximations of P RM . In some applications, the biomass, X , does not change significantly relative to S . Wenow assume that biomass is constant, i.e., X = X c . Given this assumption, (5) takesthe form T i ≤ µ max S i X c i K + S i , i ∈ N . (11)This is a convex constraint because the right hand side, a Monod function, is concave.11 heorem 3. (11) is equivalent to the hyperbolic constraint ˆ S i + ˆ T i ≤ (cid:16) ˆ S i − ˆ T i (cid:17) (cid:16) µ max KX c i + ˆ S i − ˆ T i (cid:17) (12a)0 ≤ ˆ S i − ˆ T i , (12b) where ˆ S i = µ max S i X c i and ˆ T i = KT i , i ∈ N . As with Theorem 2, this can be shown by simplifying. We refer to P R with (12) in placeof (5) and without the biomass balance, (4d), as P RM X . In standard SOC form, (12a) iswritten (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ˆ S i ˆ T i µ max KX c i (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ≤ µ max KX c i + ˆ S i − ˆ T i . Theorem 1 does not directly apply to P RM X because X = X c . The proof may beadapted by simply disregarding constraint (4d) in P R ; we omit the details because theyare straightforward. Evaluating Assumption 1 for (11) yields1 y ∂r ( s, x ) ∂s − ∂r ( s, x ) ∂x = µ max X c Ky ( K + s ) , which is always nonnegative. Therefore, (a slightly modified version of) Theorem 1always holds for P RM X . Given the approximation X = X c , P RM X is exact. As stated above, the relaxed Monod constraint (10) does not lead to a tractableoptimization problem. In this section, instead of making a physical approximation, weconstruct an SOC approximation of the non-relaxed Monod constraint, T i = µ max S i X i K + S i , i ∈ N . (13)Note that (13) is an equality, whereas (10) is an inequality. We relax this constraintusing the concave and convex envelopes in Sections 3.1 and 3.2 of [5]; see also [18, 19].We can rewrite (13) as µ max X i = T i + KT i /S i , i ∈ N . The only nonconvexity is dueto the term T i /S i . To apply the convex envelopes, we need upper and lower bounds ofthe form S i ≤ S i ≤ S i and X i ≤ X i ≤ X i for each i ∈ N . Given such bounds, we use(13), and the fact that the right hand side is increasing, to obtain the following upperand lower bounds on T i : T i = µ max S i X i K + S i , T i = µ max S i X i K + S i . (14)We know that S ≥ X ≥
0, and therefore T ≥
0. We now derive several other bounds.
Lemma 2.
For each i ∈ N , in steady state, S i ≤ max j ∈N S in j , X i ≥ min j ∈N X in j , and X i ≤ max j ∈N X in j + yS in j . Proof.
We proceed casewise.1. Starting from (4c), because T ≥ M + L ) S + CS in ≥
0, which impliesthat ( M + L ) S + C max j ∈N S in j ≥
0. Applying (2a), we have S ≤ max j ∈N S in j ,i.e., S i ≤ max j ∈N S in j for i ∈ N .2. Starting from (4d), because T ≥ M + L ) X + CX in ≤
0, which impliesthat ( M + L ) X + C min j ∈N X in j ≤
0. As above, this leads to X ≥ min j ∈N X in j ,i.e., X i ≥ min j ∈N X in j for i ∈ N .3. We combine (4c) and (4d) to obtain ( M + L ) Z + CZ in = 0, where Z = X + yS and Z in = X in + yS in . From here, similar to the previous cases, we can show thatfor each i , min j ∈N Z in j ≤ X i + yS i ≤ max j ∈N Z in j . Because S ≥
0, this implies thatfor all i ∈ N , X i ≤ max j ∈N X in j + yS in j . (cid:3) Note that if S in and X in are variables, we can simply replace them with their upperand lower bounds in Lemma 2. We can now apply the convex envelopes of [5]. For each i ∈ N , the Monod equation is represented by µ max X i = T i + Kβ i . (15a) β i represents the term T i /S i . It is constrained by the convex envelope, given below. Theconcave overestimator is the pair of linear inequalities β i SS ≤ ST i − T S i + ST (15b) β i SS ≤ ST i − T S i + ST . (15c)The convex underestimator is given by γ i ψ i ≥ T (cid:18) T − T i T − T (cid:19) (15d)( β i − γ i )( S i − ψ i ) ≥ T (cid:18) T i − TT − T (cid:19) (15e) ψ i ≥ max (cid:26) S T − T i T − T , S i − S T i − TT − T (cid:27) (15f) ψ i ≤ min (cid:26) S T − T i T − T , S i − S T i − TT − T (cid:27) (15g) β i − γ i ≥ , γ i ≥ . (15h)The first two constraints are hyperbolic SOC like (9) and (12), and the rest are linear. ψ i and γ i are auxiliary variables.We refer to P with (15) instead of (4b) as P RME . P RME may be a very good approx-imation to P , but there are no theoretical results guaranteeing exactness.13 . Linearization of bilinear terms When d , Q , and Q in are variables, the mass flow terms in constraints (4c) and (4d)are bilinear. Bilinear constraints are nonconvex and in general difficult to optimize over.Two common ways to deal with bilinearities are convex relaxations, e.g., McCormick [44]and lift-and-project [45], and disjunctive programming [16]. Disjunctive programmingtechniques have been used in several mathematically similar problems, including chemi-cal process optimization [17] and transmission network expansion in power systems [46].Here we let the flows be discrete, in which case disjunctive programming leads to ex-act linearization. We note that relaxations are more appropriate when the flows arecontinuous instead of discrete.We denote the set of pipes J . For each pipe ij ∈ J , define the binary variable λ ij ∈ { , } . (16)This could represent the decision to build a pipe, or the decision to turn on a fixed-speedpump. If both ij ∈ J and ji ∈ J , we prohibit simultaneous flow in both directions withthe constraint λ ij + λ ji ≤ . (17)The flow rate through pipe ij ∈ J is given by Q ij = Q ij + λ ij Q ij , where the constants Q ij and Q ij are respectively the base flow and the added flow if λ ij = 1. Introduce thevariable F Sij and the bilinear constraint F Sij = λ ij Q ij S i . The flow of substrate from i to j is Q ij S i + F Sij . Because λ ij ∈ { , } , we can rewrite this as the pair of linear disjunctive[16] constraints (1 − λ ij )Γ ≥ (cid:12)(cid:12) Q ij S i − F Sij (cid:12)(cid:12) , λ ij Γ ≥ (cid:12)(cid:12) F Sij (cid:12)(cid:12) , (18a)where Γ is a large positive number. Hence F Sij = 0 when λ ij = 0, and F Sij = Q ij S i when λ ij = 1. We can similarly represent the constraint F Xij = λ ij Q ij X i as the pair of linearconstraints (1 − λ ij )Γ ≥ (cid:12)(cid:12) Q ij X i − F Xij (cid:12)(cid:12) , λ ij Γ ≥ (cid:12)(cid:12) F Xij (cid:12)(cid:12) . (18b)The diffusion of substrate in pipe ij ∈ J is (cid:0) d ij + λ ij d ij (cid:1) ( S i − S j ). Introduce thevariables G Sij and G Xij and the bilinear constraints G Sij = λ ij d ij ( S i − S j ) and G Xij = λ ij d ij ( X i − X j ). We can similarly represent these as(1 − λ ij )Γ ≥ (cid:12)(cid:12) d ij ( S i − S j ) − G Sij (cid:12)(cid:12) , λ ij Γ ≥ (cid:12)(cid:12) G Sij (cid:12)(cid:12) (19a)(1 − λ ij )Γ ≥ (cid:12)(cid:12) d ij ( X i − X j ) − G Xij (cid:12)(cid:12) , λ ij Γ ≥ (cid:12)(cid:12) G Xij (cid:12)(cid:12) . (19b)If S in i and X in i are variables, then the products Q in i S in i and Q in i X in i are bilinear aswell. We linearize them by noting that, due to flow conservation, (1), and the factthat Q out is fixed, λ fully determines Q in . Let H S, ij = λ ij Q ij S in i , H S, ij = λ ji Q ji S in i ,14 X, ij = λ ij Q ij X in i , and H X, ij = λ ji Q ji X in i . These are equivalent to(1 − λ ij )Γ ≥ (cid:12)(cid:12)(cid:12) Q ij S in i − H S, ij (cid:12)(cid:12)(cid:12) , λ ij Γ ≥ (cid:12)(cid:12)(cid:12) H S, ij (cid:12)(cid:12)(cid:12) (20a)(1 − λ ji )Γ ≥ (cid:12)(cid:12)(cid:12) Q ji S in i − H S, ij (cid:12)(cid:12)(cid:12) , λ ji Γ ≥ (cid:12)(cid:12)(cid:12) H S, ij (cid:12)(cid:12)(cid:12) (20b)(1 − λ ij )Γ ≥ (cid:12)(cid:12)(cid:12) Q ij X in i − H X, ij (cid:12)(cid:12)(cid:12) , λ ij Γ ≥ (cid:12)(cid:12)(cid:12) H X, ij (cid:12)(cid:12)(cid:12) (20c)(1 − λ ji )Γ ≥ (cid:12)(cid:12)(cid:12) Q ji X in i − H X, ij (cid:12)(cid:12)(cid:12) , λ ji Γ ≥ (cid:12)(cid:12)(cid:12) H X, ij (cid:12)(cid:12)(cid:12) . (20d)After making the appropriate substitutions, we obtain the following replacements forconstraints (4c)-(4e), which we write in scalar form for each i ∈ N . For clarity, weindicate beneath each new term the corresponding entry in the original constraint. Theflow balance, (4e), becomes Q in i + X j ∈N Q ji + λ ji Q ji | {z } Q ji = Q out i + X j ∈N Q ij + λ ij Q ij | {z } Q ij . (21a)The substrate balance, (4c), becomes1 y V ii T i + Q out i S i + X j ∈N Q ij S i + F Sij | {z } Q ij S i + d ij ( S i − S j ) + G Sij | {z } d ij ( S i − S j ) = Q out i S in i + X j ∈N (cid:0) Q ij − Q ji (cid:1) S in i + H S, ij − H S, ij | {z } Q in i S in i + X j ∈N Q ji S j + F Sji | {z } Q ji S j + d ji ( S j − S i ) + G Sji | {z } d ji ( S j − S i ) . (21b)The biomass balance, (4d) takes the same form as (21b), but with X in place of S .We now make several comments. Q out i and V ii are constant. If they were insteadbinary variables, we could linearize the resulting bilinearities using the same technique.The same binary variables, λ , appear throughout (16)-(21). We could also straightfor-wardly generalize this to integer capacities by associating multiple binary variables witheach pipe.
6. Extensions
We now consider two basic extensions to P R and its SOC representations. When P R is not an exact relaxation, it is useful to constrain T from below usingan underestimator of the kinetics. To retain tractability, the underestimator should bea convex function. There are multiple ways to do this, such as using lift-and-projectrelaxations [45] or convex envelopes [18]. 15n this section, we design simple underestimators using bounds of the form S ≤ S ≤ S , X ≤ X ≤ X , and T ≤ T ≤ T . These bounds could be from Lemma 2, which makes noassumption about the growth rate except that it is nonnegative, and therefore applies forboth Monod and Contois. Alternatively, such bounds might be implied by operationalconstraints in (4f).We first consider the Contois growth rate, the kinetics of which are increasing in both S i and X i . Let T i = µ max S i X i KX i + S i , T Si = µ max S i X i KX i + S i , T Xi = µ max S i X i KX i + S i . For each i ∈ N , we have T i − T i ≥ max ( T Si − T i S i − S i ( S i − S i ) ,T Xi − T i X i − X i ( X i − X i ) ) . (22)The first argument of the maximum is a linear interpolator with X fixed at its lowerbound, and the latter with S fixed at its lower bound.Now consider the Monod growth rate. In the exact case, an underestimator is given bythe convex envelope, (15). Under the constant biomass approximation in Section 4.2.1, X = X c , the kinetics only depend on S , and we can again use the linear interpolator.Let T and T be as in (14). For each i ∈ N , we have T i − T i ≥ T i − T i S i − S i ( S i − S i ) . (23) The steady state approximation in Section 2.2 is not appropriate if the system isundergoing a transient. Instead of setting the gradostat’s derivatives to zero, we nowreplace them with numerical approximations that are linear in the variables.Suppose that there are multiple time periods, indexed by the set T and each of length∆. We index each time-varying-quantity by ( t ), e.g., S ( t ) is the substrate concentrationat time t . The objective is to maximize biogas over all time periods: X t ∈T α t F ( T ( t )) , (24)where α ∈ (0 ,
1) is a discount factor.Let D t be a numerical approximation of the derivative at time t . For example, inthe case of Euler’s explicit method with time step ∆, D t [ S ] = ( S ( t + 1) − S ( t )) / ∆. Thesubstrate and biomass balances, (4c) and (4d), become D t [ S ] = − y V T ( t ) + ( M + L ) S ( t ) + CS in ( t ) (25a) D t [ X ] = V T ( t ) + ( M + L ) X ( t ) + CX in ( t ) (25b)16or t ∈ T . Note that there are more accurate choices for D t , e.g., Runge-Kutta schemes [47].The network, as parametrized by M , L , and C , is constant, but could easily be madetime-varying as well. The remaining constraints in P R are simply enforced for all t ∈ T .Below are three scenarios modeled by this setup. • M , L , and C depend on a single vector of binary variables, η . This corresponds todesigning the system, e.g., adding new pipes, so that its performance is optimizedfor a trajectory, e.g., a sequence of operating points in a representative day. • M ( t ), L ( t ), and C ( t ) depend on a sequence of vectors of binary variables, η ( t ), t ∈ T . This corresponds to dynamically reconfiguring the system through time,e.g., choosing which valves to open or which fixed speed pumps to run in each timeperiod. • M , L , and C are constant. This corresponds to optimizing the trajectory of thesubstrate and biomass over time, subject to other operational constraints. Thephysical decisions, S in ( t ) and X in ( t ), t ∈ T , represent schedules of substrate andbiomass inflow concentrations.It may be most natural to implement the last two scenarios via receding horizon con-trol [48]. In this case, only the decisions corresponding to the first time period areimplemented. The time horizon is then pushed back by one period, the optimizationis resolved, the ‘new’ first period’s decisions are implemented, and so on. This accom-modates uncertainty by allowing the user to update the parameters, e.g., inflows andconstraints on S in ( t ) and X in ( t ), as new information becomes available.
7. Examples
We implement P RC , P RM X , and P RME on numerical examples. We solve each op-timization using the parser CVX [49] and the solver Gurobi [9]. Table 7 summarizesthe features of each optimization model. The last column refers to whether or not thecorresponding kinetics satisfy the derivative condition in Assumption 1.Model Growth rate Class Assump. 1 P Any NLP P R Any NLP P RC Contois (MI)SOCP Usually P RM Monod NLP Usually P RM X Monod, constant X (MI)SOCP Always P RME
Monod (MI)SOCP Unlikely
Table 1: Summary of optimization models
We measure the quality of the approximations in terms of the relative differencebetween the kinetics and the variable T : E = max i | r ( S i , X i ) − T i | r ( S i , X i ) . If Theorem 1 holds, then E will be zero for P RC and/or P RM X . We expect E to usuallybe positive for P RME because exactness is never guaranteed.17 .1. Variable flows
In these examples, the flows and diffusions depend on binary variables, as described inSection 5. In this case, P RC , P RM X , and P RME are MISOCPs. MISOCPs are NP-hard,but can be solved at moderate scales.We first state the models in full. The objective in each case is to maximize biogasproduction, (3), with M = N . The constraints are as follows. • A maximum budget, P ij ∈J c ij λ ij ≤ B , where c ij is the cost to install a pipe fromtank i and j , and B is the budget. This corresponds Ω in (4f). • For each i ∈ N , the SOC growth constraint. If P RC , this is (9); if P RM X , (12); if P RME , (15). In the case of P RM X , X = X c , and the other constraints on X aredropped. • For each i ∈ N , a linear underestimator of the growth constraint. If P RC , this is(22); if P RM X , (23). • Binary and linear disjuctive constraints on the flows and diffusions. For each ij ∈J , (16), (17), (18), (19). Note that we do not include (20) because S in and X in are constant in these examples. • Flow, substrate, and biomass balances. For each i ∈ N , (21). We first consider a small, four-tank example. The growth rate parameters (Monodand Contois) are µ max = K = y = 1. There is no base network, i.e., Q = d = 0. Allpairs of tanks are candidates for new pipes in either direction, so that there are twelvebinary variables. For all ij ∈ J , Q ij = 1 and d ij = 0 .
3. The cost of each new pipeis c ij = 1, and the budget is B = 4. The tank parameters are V = diag[1 2 3 4] ⊤ , Q out = [2 1 3 2] ⊤ , S in = [1 3 1 2] ⊤ , and X in = [4 3 2 1] ⊤ . In P RM X , we set X c = X in .We set Γ = 50 in all disjunctive constraints.The results are summarized in Table 7.1.1. P RC and P RM X are both exact, as pre-dicted by Theorem 1, and both result in the same pipe additions. P RME is not exact andhas a slightly different solution.Model Time (s) E Objective New pipes P RC P RM X P RME
Table 2: Results for the four-tank system
We now change the objective from P i =1 V ii T i to P i =2 V ii T i , and leave all otherparameters the same. In this case, Theorem 1 does not apply. The results are summarizedin Table 7.1.1. The pipe additions are unchanged, but now neither P RC nor P RM X areexact. In both, T binds with its underestimator, (22) or (23), and T i = r ( S i , X i ) for i = 2 , ,
4. This indicates that the assumptions of Theorem 1 are somewhat rigid, andthat when they are violated, inexactness tends to occur locally.18odel Time (s) E Objective New pipes P RC P RM X P RME
Table 3: Results for the four-tank gradostat with modified objective n tanks We now look at a larger example to see how the MISOCPs scale, using P RC as therepresentative model.The gradostat has n tanks arranged in a wheel. The first tank is a central hub, andthe other n − i ∈ J and i ∈ J for i ∈ N \ n ∈ J , n ∈ J , and i, i + 1 ∈ J and i + 1 , i ∈ J for i ∈ N \ { , n } . If there are n tanks, then there are 4 n − i ∈ N , the volumesare V ii = i and the outflows Q out i = 1. In the hard case, they are V ii = 1 + ( i mod 6)and Q out i = 1 + ( i mod 7). The other parameters, which are the same in both cases, are S in i = i , X in i = n − i + 1, and maximum budget B = 1 . n , and the rest are the same asin the previous example.Figure 1 show the results. In the easy case, the computation time increases roughlylinearly, taking around a minute with with 60 tanks and 236 binary variables. In the hardcase, the computation time increases exponentially, and more than an hour is needed with11 tanks and 44 binary variables.The reason for the difference is that in the easy case, tanks with larger volumesare significantly better represented in the objective. This enables the solver to rapidlyeliminate solutions with many pipes added to the smaller tanks, whereas in the hardercase, the solver must search the feasible set more evenly. P RC is exact in all cases. We now implement P RC on an example with multiple time periods, as described inSection 6.2. The water flows and diffusions are constant, which makes P RC an SOCP.Our intention here is to demonstrate that this version of the problem can be solved atvery large scales. This enables us to deal with larger systems and, as in this example,to choose the time step small enough that the continuous dynamics of the gradostat arewell-represented.There are four tanks, all with unit volume and growth rate parameters. The inflowvector is Q in = [2 1 1 1] ⊤ . The flows between tanks are: Q = 1, Q = 2, Q = 1, Q = 1, and the diffusion is d = 0 . Q .There are τ = 1000 time periods of length ∆ = 1. We approximate the derivative withEuler’s method. The inflow substrate concentrations, shown in the top plot of Figure 2,are S in1 ( t ) = 1 + sin(4 πt/τ ), S in2 ( t ) = 0, S in3 ( t ) = (cid:26) / , τ / < t ≤ τ / , otherwise ,
20 40 60 80 100 120 140 160 180 200 220
Number of binary variables S e c ond s Figure 1: Computation time versus number of binary variables for P RC , an MISOCP, on a wheelgradostat. The hard case is shown on top and the easy case on bottom. Note that the top plot has alogarithmic y -axis. and S in4 ( t ) = 1 + cos(4 πt/τ ) for t ∈ T .The objective is to maximize the cumulative biogas production over all time periods,(24). The constraints are listed below. • At each time t ∈ T , the total biomass added must satisfy Q in ⊤ X ( t ) ≤ • For each i ∈ N and t ∈ T , the SOC Contois growth constraint, (9). • At each time t ∈ T , the dynamic substrate and biomass balances, (25). Theseconstraints couple the variables in consecutive time periods. • As boundary conditions we require that S (1) = S ( τ + 1) and X (1) = X ( τ + 1).There are roughly sixteen thousand variables in this problem. The computation timewas 281 seconds, of which only 1.49 were taken by the solver and the rest by the parser.The optimal objective was 1140.18 units of biogas mass.Figure 2 shows S in ( t ), which is as specified above, and the optimal X in ( t ). Figure 3shows the optimal S ( t ), X ( t ), and T ( t ). 20
00 200 300 400 500 600 700 800 900 100000.511.52 S i n ( t ) Tank 1Tank 2Tank 3Tank 4
100 200 300 400 500 600 700 800 900 1000
Time step X i n ( t ) Figure 2: S in ( t ) and X in ( t ) as functions of time X in ( t ) is zero in tanks 2 and 3, except after S in3 ( t ) jumps up to 1 / X in ( t ) and the curves in Figure 3to capture the oscillations caused by these transitions.While S in3 ( t ) = 1 /
2, the substrate, biomass, and production in all tanks except thefirst jump, although slightly in tanks 2 and 4. This is because the additional substrateinjected into tank 3 reaches tanks 2 and 4 through the flows. Tank 1 only receivessubstrate and biomass from tank 2 through a small diffusive coupling, and for this reasonis not noticeably affected by the change in S in3 ( t ).The solution was exact in all times periods. However, we observed that when therewere fewer time periods and a larger time step, the solution was not always exact. Wealso remark that trajectories produced by P RC do not necessarily lead to better stabilityor disturbance rejection. Such control objectives could be incorporated through a track-ing objective and a receding horizon implementation. Characterizing the exactness ofdynamic optimizations with different objectives is a topic of future work.
8. Conclusions
We have formulated SOCPs for optimizing the gradostat with Contois or Monodgrowth rates. The SOCPs are convex relaxations, which we proved are exact under21
00 200 300 400 500 600 700 800 900 10000.20.40.60.811.21.4 S ( t )
100 200 300 400 500 600 700 800 900 10000.511.5 X ( t )
100 200 300 400 500 600 700 800 900 1000
Time step T ( t ) Figure 3: S ( t ) and X ( t ) as functions of time simple conditions. We also gave linear underestimators, which are useful when the re-laxations are not exact, and a dynamic extension in which the derivatives are replacedwith numerical approximations instead of set to zero.There are many directions for future work. More physical features could be in-corporated, including continuous variable flows via convex relaxations, recirculation ofbiomass, multi-reactions with several substrates, and other growth rates such as Teissierand Haldane. The constraints of the SOCPs could be used for a number of other pur-poses such as state estimation, setpoint tracking, and receding horizon control. As inthis paper, the KKT conditions could be used to identify conditions under which each ofthese modifications is exact. Acknowledgements
We thank Professor Denis Dochain for helpful discussion.
References [1] Y. Shen, J. Linville, M. Urgun-Demirtas, M. Mintz, S. Snyder, An overview of biogas productionand utilization at full-scale wastewater treatment plants (WWTPs) in the United States: challenges nd opportunities towards energy-neutral WWTPs, Renewable and Sustainable Energy Reviews 50(2015) 346–362.[2] J. Monod, The growth of bacterial cultures, Annual review of microbiology 3 (1) (1949) 371–394.[3] D. Contois, Kinetics of bacterial growth: relationship between population density and specificgrowth rate of continuous cultures, Microbiology 21 (1) (1959) 40–50.[4] M. Lobo, L. Vandenberghe, S. Boyd, H. Lebret, Applications of second-order cone programming,Linear Algebra and its Applications 284 (1998) 193–228.[5] M. Tawarmalani, N. Sahinidis, Semidefinite relaxations of fractional programs via novel convexifi-cation techniques, Journal of Global Optimization 20 (2) (2001) 133–154.[6] S. Drewes, Mixed integer second order cone programming, Ph.D. thesis, Technischen Universit¨atDarmstadt, Department of Mathematics (2009).[7] A. Atamt¨urk, V. Narayanan, Conic mixed-integer rounding cuts, Mathematical Programming 122(2010) 1–20.[8] P. Belotti, J. G´oez, I. P´olik, T. Ralphs, T. Terlaky, A complete characterization of disjunctive coniccuts for mixed integer second order cone optimization, Discrete optimization 24 (2017) 3–31.[9] Gurobi Optimization, LLC, Gurobi optimizer reference manual (2021).URL [10] H. Smith, P. Waltman, The theory of the chemostat: dynamics of microbial competition, Vol. 13,Cambridge University Press, 1995.[11] J. Harmand, C. Lobry, A. Rapaport, T. Sari, The chemostat: Mathematical theory of microorganismcultures, John Wiley & Sons, 2017.[12] R. Lovitt, J. Wimpenny, The gradostat: a bidirectional compound chemostat and its applicationin microbiological research, Microbiology 127 (2) (1981) 261–268.[13] W. J¨ager, J.-H. So, B. Tang, P. Waltman, Competition in the gradostat, Journal of MathematicalBiology 25 (1) (1987) 23–42.[14] T. Bayen, A. Rapaport, M. Sebbah, Minimal time control of the two tanks gradostat model undera cascade input constraint, SIAM Journal on Control and Optimization 52 (4) (2014) 2568–2594.[15] L. Biegler, I. Grossmann, A. Westerberg, Systematic methods of chemical process design, PrenticeHall, 1997.[16] I. Grossmann, Review of nonlinear mixed-integer and disjunctive programming techniques, Opti-mization and Engineering 3 (3) (2002) 227–252.[17] S. Lee, I. Grossmann, Global optimization of nonlinear generalized disjunctive programming withbilinear equality constraints: applications to process networks, Computers & Chemical Engineering27 (11) (2003) 1557–1575.[18] M. Tawarmalani, N. Sahinidis, Convex extensions and envelopes of lower semi-continuous functions,Mathematical Programming 93 (2) (2002) 247–263.[19] M. Tawarmalani, N. Sahinidis, Convexification and global optimization in continuous and mixed-integer nonlinear programming: theory, algorithms, software, and applications, Vol. 65, SpringerScience & Business Media, 2013.[20] G. Froment, K. Bischoff, J. De Wilde, Chemical reactor analysis and design, 3rd Edition, Wiley,2010.[21] R. Aris, The Optimal Design of Chemical Reactors: A Study in Dynamic Programming, AcademicPress, 1962.[22] K. Bischoff, Optimal continuous fermentation reactor design, The Canadian Journal of ChemicalEngineering 44 (5) (1966) 281–284.[23] K. C. A. Luyben, J. Tramper, Optimal design for continuous stirred tank reactors in series usingMichaelis–Menten kinetics, Biotechnology and Bioengineering 24 (5) (1982) 1217–1220.[24] G. Hill, C. Robinson, Minimum tank volumes for CFST bioreactors in series, The Canadian Journalof Chemical Engineering 67 (5) (1989) 818–824.[25] C. de Gooijer, W. Bakker, H. Beeftink, J. Tramper, Bioreactors in series: an overview of designprocedures and practical applications, Enzyme and Microbial Technology 18 (3) (1996) 202–219.[26] J. Harmand, A. Rapaport, A. Trofino, Optimal design of interconnected bioreactors: New results,AIChE Journal 49 (6) (2003) 1433–1450.[27] J. Harmand, A. Rapaport, A. Dram´e, Optimal design of two interconnected enzymatic bioreactors,Journal of Process Control 14 (7) (2004) 785–794.[28] J. Harmand, D. Dochain, The optimal design of two interconnected (bio) chemical reactors revisited,Computers & chemical engineering 30 (1) (2005) 70–82.[29] J. Zambrano, B. Carlsson, Optimizing zone volumes in bioreactors described by Monod and Contoisgrowth kinetics, in: Proceeding of the IWA World Water Congress & Exhibition, Lisbon, Portugal, doi:https://doi.org/10.1016/j.ifacol.2018.09.236 .[34] C. Robles-Rodriguez, A. Ben-Ayed, J. Bernier, V. Rocher, D. Dochain, Management of an integratednetwork of wastewater treatment plants for improving water quality in a river basin, in: IFACSymposia Series, Vol. 52, 2019, p. 358.[35] G. Bastin, D. Dochain, On-line estimation and adaptive control of bioreactors, Vol. 1, Elsevier,2013.[36] S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.[37] J. Jacquez, C. Simon, Qualitative theory of compartmental systems, SIAM Review 35 (1) (1993)43–79.[38] A. Berman, R. Plemmons, Nonnegative matrices in the mathematical sciences, SIAM, 1994.[39] A. Rapaport, Lecture notes on mathematical models of interconnected chemostats (Jul. 2019).URL https://hal.archives-ouvertes.fr/cel-02182956 [40] J. Schimel, M. Weintraub, The implications of exoenzyme activity on microbial carbon and nitrogenlimitation in soil: a theoretical model, Soil Biology and Biochemistry 35 (4) (2003) 549–563.[41] W. Parton, J. Stewart, C. Cole, Dynamics of C, N, P and S in grassland soils: a model, Biogeo-chemistry 5 (1) (1988) 109–131.[42] D. Moorhead, R. Sinsabaugh, A theoretical model of litter decay and microbial interaction, Eco-logical Monographs 76 (2) (2006) 151–174.[43] D. Moorhead, G. Lashermes, R. Sinsabaugh, A theoretical model of C-and N-acquiring exoenzymeactivities, which balances microbial demands during decomposition, Soil Biology and Biochemistry53 (2012) 133–141.[44] G. McCormick, Computability of global solutions to factorable nonconvex programs: Part I—Convex underestimating problems, Mathematical programming 10 (1) (1976) 147–175.[45] H. Sherali, A. Alameddine, A new reformulation-linearization technique for bilinear programmingproblems, Journal of global optimization 2 (1992) 379–410.[46] L. Bahiense, G. Oliveira, M. Pereira, A mixed integer disjunctive model for transmission networkexpansion, IEEE Transactions on Power Systems 16 (2001) 560–565.[47] J. Betts, Survey of numerical methods for trajectory optimization, Journal of Guidance, Control,and Dynamics 21 (2) (1998) 193–207.[48] J. Mattingley, Y. Wang, S. Boyd, Receding horizon control, IEEE Control Systems Magazine 31 (3)(2011) 52–65.[49] M. Grant, S. Boyd, CVX: Matlab software for disciplined convex programming, version 2.1, http://cvxr.com/cvx (Mar. 2014).(Mar. 2014).