Optimization of Fractional Order Dynamic Chemical Processing Systems

Article pubs.acs.org/IECR

Optimization of Fractional Order Dynamic Chemical Processing Systems Antonio Flores-Tlacuahuac* and Lorenz T. Biegler Department of Chemical Engineering Carnegie-Mellon University, 5000 Forbes Avenue, Pittsburgh, Pennsylvania, 15213, United States ABSTRACT: In this work we address the dynamic simulation and optimization of chemical processing systems modeled in terms of f ractional order dif ferential equations. While fractional derivatives were first proposed by Liouville in 1832 [Samko et al. Fractional Integrals and Derivatives Theory and Applications; Gordon and Breach: New York, 1993; Oldham and Spanier. The f ractional Calculus; Academic Press: New York, 1974], fractional differential equation (FDE) models have been only recently been explored. These have been proposed for a wide range of applications that include systems with nonlocal diffusion phenomena and geometries with fractal dimensions. FDE models have been shown to have advantages over traditional integer order models, as they often avoid scale dependent model parameters. For medium or large scale applications of FDEs normally no analytical solutions are available, and therefore, approximated numerical solutions ought to be sought. Moreover, little work has been done to solve fractional order differential equations numerically; most of the existing numerical methods are intended for small scale systems. In this work, we propose a new numerical method, based on Gaussian quadrature on finite elements, to address larger-scale FDEs and extend them to dynamic optimization. To gain a better appreciation about the performance of the new algorithm, we compare its response to recently proposed predictor−corrector methods. We also develop a proposed method to deal with dynamic optimization fractional order systems, an open research problem that has not received wide consideration. We test the performance of the algorithms by deploying three systems embedded with different fractional derivative behavior, from a simple linear dynamic system to a dynamic multiple steady-states bioreactor with various levels of imperfect mixing. The results indicate better numerical properties of the fractional order Gauss quadrature algorithm both for the dynamic simulation and optimal control issues.

1. INTRODUCTION Mathematical models of macrochemical processing systems are ubiquitous in Chemical Engineering. On the other hand, such models need to be extended to include nanoscale behavior and to take into account the interactions among multiscale phenomena. In particular, it has been reported that certain kinds of distributed or dynamic natural phenomena, such as nonlocal diffusion and complex geometries modeled in fractal dimensions, cannot be properly modeled in terms of traditional integer order differential equations.1−5 These complex, distributed behaviors have been encountered in viscoelastic systems,6 diffusion,7 transport phenomena,8 control,9 market dynamics,10 fractional kinetics,11 pharmacokinetics,12,13 DNA code analysis,14 advection-dispersion systems in earth science,15 tumor growth,16 and even in the social sciences.17 Modeling successes for these phenomena have recently been reported using fractional derivatives for both ordinary and partial differential equations. The role of fractional derivatives in modeling complex physical phenomena has been explored in several recent studies. Fractional derivatives are often used to model nonlocal transport and this behavior frequently arises in complex microstructural systems that can be modeled with fractal geometries. However, the modeling of this behavior differs from application to application, depending on whether the modeling framework is based on particle approximations and stochastic analysis of non-Gaussian distributions, or aggregation procedures used to model microstructures with fractal dimensions. Nevertheless, fractional differential operators arise for both of © 2013 American Chemical Society

these phenomena in the following way. For both of these phenomena, fractional differential equations (FDEs) can be modeled in the Laplace domain through continued fraction expansions and the aggregation of these expansions (under assumptions of constant microstructural parameters) leads to fractional exponents in the Laplace domain. In particular, we consider three important studies that provide detailed and complex, physical interpretations of fractional derivatives. Schumer et al.15 developed FDE models to describe nonlocal dispersion behavior in time and space. Derived from a stochastic diffusion framework, fractional derivatives can be used to describe non-Fickian dispersion with memory effects. These nonlocal effects have are due to particle distributions with nonGaussian spreads, e.g., random walk behavior in time and through particle hops with delays, leading to non-Fickian diffusion, with spatial fractional derivatives between one and two. These phenomena have been observed in particle flows and flow through porous media. In particular, this particle distribution analysis can predict superdiffusion due to particle jump or hop behaviors, leading to faster spreading than Fickian diffusion. On the other hand, subdiffusion can be modeled for particle flows that manifest dead zones with slower distribution Special Issue: David Himmelblau and Gary Powers Memorial Received: Revised: Accepted: Published: 5110

April 25, 2013 December 27, 2013 December 31, 2013 December 31, 2013 dx.doi.org/10.1021/ie401317r | Ind. Eng. Chem. Res. 2014, 53, 5110−5127

Industrial & Engineering Chemistry Research

Article

to deal with fractional systems and to solve the resulting equations numerically.30−35 These well-known classical procedures36 lead to two-point boundary problems, which are not easy to extend beyond small scale systems. For medium or large scale systems, a number of optimization algorithms for IODES37 can also be considered for FDEs. In this work we transform the dynamic mathematical models of some well-known systems (flash units, stirred tank reactors) from integer to fractional order systems, and therefore equip them with nonlocal diffusion and incomplete mixing phenomena, as observed in our analysis of their open-loop dynamic behavior. We apply a full discretization scheme of the FDE based on the predictor-corrector method of Diethelm.22,38,24 Here we deploy a numerical predictor-corrector method24 and formulate a new numerical integration method based on approximating a Volterra equation naturally occurring in the solution of FDEs. In addition, we develop a new discretization approach based on Gauss quadrature and related to orthogonal collocation on finite elements.37 This leads to a new contribution for dynamic optimization of nonlinear FDEs proposed for chemical processing systems. Using direct transcription for the fully discretized Volterra equations, we formulate an efficient nonlinear programming (NLP) strategy for the dynamic optimization of FDE systems. Although we do not develop or interpret FDE models in the context of chemical processing systems, we hope this work will aid in its development by developing efficient numerical methods for the simulation and optimization of FDE-based systems. In the next section, we present key definitions and background on fractional calculus and the formulation of FDEs. Section 3 then develops the fully discretized predictorcorrector method for FDEs as well as a discretization approach based on Gaussian quadrature. The extension of these approaches to dynamic optimization problems is also discussed in this section. Section 4 provides three case studies, a simple linear system, a flash unit, and a bioreactor. These three studies demonstrate the nature of the FDE solution as well as the performance of our numerical methods for simulation cases. In addition, several optimal dynamic state transitions are considered to demonstrate the direct extension of these formulations to numerical solution of optimal control problems. Section 5 concludes the paper and discusses areas for future work.

spreads than Fickian diffusion. As a result, fractional derivatives in space can model high-velocity contrast in convectiondiffusion systems, while fractional derivatives in time can model nonideal residence times. Moreover, FDEs were also derived by Ionescu,13 who described branched flow behavior in self-similar networks with geometries described by fractal dimensions. Applied to transport and storage of glycogen molecules, FDEs were developed by modeling (infinite) ladder network expansions that directly lead to determination of the fractional derivative. Similarly, vascular microstructure (represented by Koch trees) leads to an infinite expansion of branched fractal structures. In both cases, using the method of continued fractions and collapsing then in the Laplace domain leads to fractional derivatives between zero and one; the derivative value can be determined directly from the frequency and lengths of the branched chains. Du et al.18 consider the analysis, solution, and approximation of nonlocal diusion problems in spatial and temporal dimensions. They develop variational principles for nonlocal vector calculus and relate analogs between nonlocal and classical diffusion models. In particular, they show that fractional Laplacian and derivative models are special cases of nonlocal diusion. As in ref 15, a probabilistic perspective shows that the nonlocal spatial operator corresponds to the infintesimal generator for a jump process. Moreover, they consider the imposition of volume constraints, which leads to a well-posed problem for fractional derivatives between zero and one. They also develop an analytical framework that allows finitedimensional approximations using discontinuous and continuous Galerkin methods, along with an analysis of error conditions. Results of this work have applications of fractional derivative equations that range from continuum mechanics to graph theory. On the other hand, we note that determining the exact value of the fractional derivative requires detailed (and often idealized) information on microstructure or distributional behavior. Since the resulting FDE model may not be easy to realize based on this information alone, the value of the fractional derivative is often used as a tuning parameter in many studies, where the model is based on experimental data. The mathematical concept of fractional derivatives, called f ractional calculus, dates back to a paper by Liouville in 1832.60 Since then, theoretical properties and mainly analytical solution methods have been advanced in this area. In contrast, numerical algorithms to solve fractional order differential equations (FDEs)19 have lagged the development of corresponding methods for integer order differential equations (IODE), which include the solution of stiff and nonstiff systems,20 along with similar advances for the numerical integration of differential and algebraic equations.21 In recent years, there have been a few attempts to extend IODE numerical methods (Euler, Runge−Kutta, Predictor-Corrector) to FDEs,22−24 but currently only small scale FDE systems have been addressed with existing numerical integration methods. Similarly, dynamic optimization of systems modeled by IODEs have witnessed enormous advances that have led to numerous process applications including grade transitions in polymerization reactors,25,26 optimal interactions between design and control,27,28 and the nonlinear model predictive control of large scale systems.29 On the other hand, studies on the optimization of FDE systems tend to extend classical Hamiltonian-based optimality conditions for dynamic systems

2. FDES: THEORETICAL BACKGROUND AND PROPERTIES We first provide an overview of some common properties of FDEs that will set the stage for numerical algorithms of such systems. FDEs can be represented either in the Riemann− Liouville or Caputo forms depending on the way initial conditions are specified. More discussion on the background and theoretical properties of FDEs can be found in refs 39 and 1−4. The traditional Riemann−Liouville definition of the fractional derivative is given by Γ(δ) dα δ x δ−α αx = Γ(δ − α + 1) dx

(1)

where δ is a noninteger exponent, α is the fractional order of the derivative term, and Γ stands for the gamma function. More generally, fractional order derivatives can be developed by applying the Cauchy formula for repeated integration: 5111

dx.doi.org/10.1021/ie401317r | Ind. Eng. Chem. Res. 2014, 53, 5110−5127

Industrial & Engineering Chemistry Research t

∫a ∫a =

s1

···

∫a

sm − 1

1 (m − 1)!

Article

This forms the basis of the numerical methods developed in the next section.

f (sm) dsm....ds1

∫a

t

(t − s)m − 1f (s) ds

3. APPROACHES FOR THE NUMERICAL SOLUTION OF FDES Only a few algorithms have been proposed for numerical integration of nonlinear FDEs. In this section we describe the predictor-corrector method proposed in refs 22 and 24 for solution of nonstiff or mildly stiff FDEs. We also provide the derivation of a new numerical algorithm for numerical solution of FDEs based on Gaussian quadrature on finite elements. Through this development we demonstrate the performance of these methods with numerical experiments and results. These help to motivate the analysis of theoretical properties (stability, convergence, approximation errors, etc), which will be left for future study. 3.1. Predictor-Corrector Method. We present a predictorcorrector method for solving fractional order differential equations (FDEs) by first recalling the corresponding derivation for IODEs. Integer Order Differential Equations. Given the following first-order initial value problem

(2)

which yields Jma f(t), the mth antiderivative (see ref 40, p 193). The right-hand side of the above equation can be cast as Jaα f (t ) =

1 Γ(α)

∫a

t

(t − s)α− 1f (s) ds

(3)

where we normally take the initial time a = 0. This is known as the Riemann−Liouville integral, which along with Dαf(t) = J−αf(t) and J−α−mJmf(t) = J−αf(t), allows us to use the following definition of a fractional derivative as def

D α f (t ) = J m − α D m f (t )

(4)

Here m = ⌈α⌉ is the smallest integer greater than α, and D = dm/dxm is the normal integer derivative operator. In the above definition, we also assume that f(t) is continuous and m times differentiable. To derive numerical methods for FDEs, we first consider a single FDE of order α as follows: m

Dα x(t ) = f (x(t ), t )

(5)

Dx(t ) = f (x(t ), t )

(11)

x(0) = y0

(12)

41

Using the Caputo representation, the right-hand side of (5) is given by Dα x(t ) = J m − α x(m)(t )

we first consider the time interval [0, T] which has been divided into an equally spaced grid of points. To compute the solution at a given point tn+1, given x(tn), we integrate 11 between these limits:

(6)

where x(m)(t) represents the mth derivative of x with respect to t (i.e., x(m)(t) = Dmx(t)). Before we consider algorithms for solution of FDEs, some discussion on specifying the initial conditions is needed. There has been some controversy about the most convenient way of specifying initial conditions in fractional systems.42 For the FDE in Riemann−Liouville form (5), the linear operator Dα has an m-dimensional kernel. To find the unique solution to the FDE, m initial conditions need to be specified. From the theory of fractional systems,43 the initial conditions are given by dα − k x(t )|t = 0 + = bk , dt α − k

k = 0, ... m − 1

x(tn + 1) = x(tn) +

k = 0, 1, ..., m − 1

∫t

m−1

∑ x0k k=0

tk 1 + Γ(α) k!

∫0

t

(13)

tn + 1

f (x(t ), t ) dt

n

(7)

≈

⎛ tn + 1 − tn ⎞ ⎜ ⎟[f (x(t ), t ) + f (x(t n n n + 1), tn + 1)] ⎝ ⎠ 2

(14)

where x(tn + 1) = x(tn) +

h [f (x(tn), tn) + f (x(tn + 1), tn + 1)] 2

(15)

and h = tn+1 − tn is the integration step size. If the function f(x(t),t) is nonlinear, then x(tn+1) must be determined by solving a system of nonlinear equations. To reduce this computational load, we approximate the value of x(tn+1) by xP(tn+1) and use the following equation:

(8)

(9)

x(tn + 1) = x(tn) +

In this work we apply the Caputo interpretation of fractional differential equations to specify the initial conditions. Moreover, from (8) as well as (3), we recognize that the solution to this FDE now satisfies the following Volterra integral equation: x(t ) =

f (x(t ), t ) dt

The value of the integral in the right-hand side of the above equation can be approximated by a number of implicit quadrature formulas. Among the simplest of these is the trapezoidal quadrature:

where Pm−1[x] is an m-order Taylor polynomial of x(t) expanded around t = 0, and the initial conditions are the ones normally used for integer order differential equations: x(k)(0) = x0(k) ,

tn + 1

n

where bk are numerical values of certain fractional derivatives. Often, the initial values for FDEs are hard to determine and their physical meaning also might not be clear. To provide clearer definition, Caputo41 suggested the specification of classical integer order derivatives for the initial conditions by casting the fractional system as Dα (x − Pm − 1[x])(t ) = f (x(t ), t )

∫t

h [f (x(tn), tn) + f (x P(tn + 1), tn + 1)] 2 (16)

P

where x (tn+1) can be computed using a rectangular quadrature approximation x P(tn + 1) = x(tn) + hf (x(tn), tn)

(t − s)α − 1f (x(s), s) ds

(17)

This leads to the explicit improved Euler (or Heun’s) method, which is suitable for nonstiff or mildly stiff dynamic systems.

(10) 5112



Article

Fractional Order Differential Equations. For the numerical solution of the single FDE: α

D x(t ) = f (x(t ), t )

(18)

x k(0) = x0k ,

(19)

bj , n + 1 =

m−1

k = 1, ..., m

tn + 1

(tn + 1 − s)α − 1g (s) ds =

∫0

tn + 1

x P(tn + 1) =

∫0

(tn + 1 − s)

m−1

x(t ) =

+

∑ aj ,n+ 1g(t j) j=0

⎧ nα + 1 − (n − α)(n + 1)α , if j = 0 ⎪ α+1 ⎪ + (n − j)α + 1 = ⎨(n − j + 2) if 1 ≤ j ≤ n ⎪ − 2(n − j + 1)α + 1 , ⎪ if j = n + 1 ⎩1,

∫t

t

k=0

+

∫t

t

∫0

tn ,0

(t − s)α − 1f (x(s), s) ds

(t − s)α − 1f (x(s), s) ds

n ,0

t

(t − s)α =− f (x(s), s) α t

+

∫t

t

n ,0

n ,0

(28)

(t − s)α α

((∂f /∂x)T x ̇ + ∂f /∂s) ds

(29)

1 (t − s)α− 1f (x(s), s) Γ(α)

(30)

where

tnk+ 1 (k) hα x0 + f (x(tn + 1), tn + 1) Γ(α + 2) k!

hα Γ(α + 2)

1 Γ(α)

1 Γ(α)

n ,0

Therefore, the integral in the right-hand side of (10) can be approximated by

∑

k!

+

(t − s)α − 1f (x(s), s) ds

(22)

m−1

k

where tn,0 is the left boundary of finite element n ∈ {1, ..., N} and t ∈ [tn,0, tn+1,0]. The first integral can be evaluated using kpoint Gaussian quadrature over finite elements, as the integrand is well-defined for all values of α. Applying a quadrature formula to the second integral is more difficult; if a “closed” formula is used, then (t − s)α−1 cannot be evaluated at s = t when α < 1. As a result, we apply integration by parts to the second term so that

n+1

where

x(tn + 1) =

∑ x0k t k=0

(21)

aj , n + 1

j=0

Since this method is based on a second order quadrature formula, we expect the resulting predictor-corrector method to have a global error of O(h2). 3.2. Algorithm Based on Collocation on Finite Elements. To develop a higher order, implicit approach that can be extended to stiff systems, we recall the Volterra integral equation (10) and rewrite this equation as

(tn + 1 − s)α − 1gn + 1(s) ds

hα gn + 1(s) ds = α(α + 1)

∑ bj ,n+ 1f (x(t j), t j) (27)

where gn+1(s) is a piecewise linear interpolant of g(s) evaluated at the nodes tj,j = 0,1, ..., n + 1. The integral in the right-hand side of the above equation reads as follows: α−1

n

tnk+ 1 (k) 1 x0 + Γ(α) k!

∑ k=0

(20)

tn + 1

(26)

Therefore, the predictor equation (PE) reads as follows:

we consider and apply the predictor-corrector method developed in refs 22 and 44. This uses an equation similar to (13) based on (10). However, the integrand in (10) is a product of functions, so a different treatment is required to approximate its value compared to numerical integration of single functions. Here numerical evaluation of the integral of the product of functions is done with a product trapezoidal quadrature formula:

∫0

hα ((n + 1 − j)α − (n − j)α α

ẋ =

leaving

n

∑ aj ,n+ 1f (x(t j), t j)

∫t

(23)

j=0

t

(t − s)α − 1f (x(s), s) ds

n ,0

As in (14) x(tn+1) appears on both sides of the above equation and we again use an approximate value of the solution x(tn+1). Hence, we obtain the corrector equation (CE): m−1

x(tn + 1) =

∑ k=0

=

tn + 1

m−1

∑ aj ,n+ 1f (x(t j), t j)

x(t ) =

+

(t − s)α α Γ(α)

T

(∂f /∂x) f (x(s), s) + ∂f /∂s) ds

(31)

k

k!

+

1 Γ(α)

∫0

tn ,0

(t − s)α − 1f (x(s), s) ds

α 1 ⎡ (t − tn ,0) ⎢ f (x(tn ,0), tn ,0) α Γ(α) ⎣

(t − s)α ((t − s)α − 1(∂f /∂x)T f (x(s), s) ( ) α Γ α n ,0 ⎤ + ∂f /∂s) ds ⎥ ⎦ (32) +

∑ bj ,n+ 1g(t j) j=0

∑ x0k t k=0

(24)

n

(tn + 1 − u)α − 1g (u) du =

t

The resulting equation is given by

n j=0

∫t

n ,0

((t − s)

To compute the predictor solution xP(tn+1), we once again use the product rectangle quadrature formula to approximate the integral term on the right-hand side of (10):

∫0

α

f (x(tn ,0), tn ,0) +

α−1

tnk+ 1 (k) hα x0 + f (x P(tn + 1), tn + 1) Γ(α + 2) k!

hα + Γ(α + 2)

(t − tn ,0)α

(25)

where 5113

∫t

t



Article

where the integrand is now well-defined for α > 1/2. Notice that for smaller values of α, additional integrations by parts would be required. The second integral in (32) is now evaluated using a trapezoidal method applied at k Gauss points tnj in the interval n, thus leading to m−1

x(tnj) =

∑ x0k t

k

k!

k=0

+

1 Γ(α)

n−1

approach all the system states, input, and output variables are discretized at the same level as the state variables through piecewise polynomial profiles. Because today this is well-known material, we did not include details on how this discretization is approached. More information on how to address the discretization issue can be found in ref 45.

k

∑ ∑ (tnj − tn′j′)α − 1f

4. EXAMPLES In this section we apply the two numerical methods described in section 3 to the dynamic simulation and optimization of chemical processing systems. We begin with a simple linear example, which has an analytical solution that can be used to check the accuracy of the developed numerical methods. We then consider to two process systems featuring nonlinear behavior and also modified to FDEs to consider nonlocal diffusion and nonideal residence times. Although dynamic optimization techniques have been widely applied to a large number of complex systems, the application of optimal control to such fractional models has been rather limited with just a a few small scale published examples.34 Initial efforts on the application of optimal control to fractional systems have concentrated on extending the classical optimal control theory based on Pontryagin maximium principle to fractional systems.30−33,35 Even though some interesting open-loop dynamic simulation results have recently been published,46−51 to our best knowledge, no results on the optimal control of fractional order systems have been published for chemical processing systems. Our aim in the present work is to deal with the optimal control of small scale fractional order systems featuring embedded nonlinearities, such that the advantages of the two proposed optimal control algorithms can be compared and appreciated, and eventually extended to cope with the dynamic optimization of larger scale complex systems. All of the examples below were modeled in AMPL52 or GAMS,53 using IPOPT54 as the NLP solver. Even with the finest grid and smallest stepsizes, all of the cases presented could be solved reasonably quickly, within 1 CPU min (2.66 GHz Intel 2-Core processor). 4.1. Fundamental Fractional Linear Differential Equation. As the first case study for the numerical solution of fractional order differential systems, we use the so-called fundamental linear fractional order differential equation:31

n′= 1 j′= 1

α ⎡ 1 ⎢ (tnj − tn ,0) (x(tn′j′), tn′j′) + f (x(tn ,0), tn ,0) Γ(α) ⎢⎣ α j

+

∑ j′= 1

(tnj′ − tnj′− 1) 2α Γ(α)

((∂f /∂x)T {(tnj − tnj′− 1)2α − 1f

(x(tnj′− 1), tnj′− 1) + (tnj − tnj′)2α − 1f (x(tnj′), tnj′)} + {(tnj − tnj′− 1)α ∂f (x(tnj′− 1), tnj′− 1)/∂t + (tnj − tnj′)2α − 1f (x(tnj′), tnj′)} + {(tnj − tnj′− 1)α ⎤ ∂f (x(tnj′− 1), tnj′− 1)/∂t + (tnj − tnj′) ∂f (x(tnj′), tnj′)/∂t })⎥ ⎥⎦ α

(33)

Moreover, since this method is based on a k-point Gaussian quadrature formula, we expect the resulting collocation method to have a global error of O(h2k). As a result, for an equivalent accuracy, far fewer time steps (finite elements) are required than with the predictor-corrector method. 3.3. Extension to Dynamic Optimization via Direct Transcription. We now consider the solution of dynamic optimization problems with FDEs of the form: min Φ(x(t f ))

(34a)

s.t. Dαi x = fi (x(t ), u(t ), p , t ), xiki(0) = xik,0i ,

i = 1, ..., nx

(34b)

i = 1, ... nx , k = 0, ..., ⌈αi ⌉ − 1

(34c)

Here x(t) ∈  are the state variables and u(t) ∈  and p ∈ np are control profiles and time-independent optimization variables, respectively. To solve this problem, we apply either the predict-evaluate-correct-evaluate equations (PECE) or the Gaussian quadrature discretization described earlier in this section and formulate the following nonlinear programming (NLP) problem. After discretizing the state and control profiles and substituting the PECEs (24) and (27) or the quadrature eqs 33 into section 4.1, we obtain: nx

nu

Dα x = −ax(t ) + bu(t )

where x is the state variable, t is the independent variable, u is a control variable, α is the order of the equation, and a, b are constants with a = b = 1. Fractional System Dynamic Simulation. Figure 1 displays the open-loop dynamic response of the addressed system using an arbitrary constant value of the control action u = −0.3. As the fractional order of the model becomes smaller we see that the system exhibits a slower time response. The results shown in Figure 1a were computed using our own implementation of the PECE algorithm described previously22,38,24 and double checked first by converting the time domain model into the Laplace domain and then numerically inverting it55 to get the time domain behavior. Results were also checked against a recent PECE implementation described in ref 44. In addition, results from the proposed Gaussian quadrature discretization approach are shown in Figure 1b. Although not shown, the dynamic simulation results obtained using the PECE and Gaussian quadrature algorithms match well with the Laplace solution and the PECE implementation in ref 44.

min Φ(xN ) s.t.

discretized FDEs (24), (27), or (33)

(36)

(35)

Problem 35 allows for a number of formulation options. In particular, the fractional order can be different for each differential equation, and corresponding initial conditions can be specified as parameters for optimization. For this formulation, we note that the state profiles are described by a continuous discretization over time while control profiles are not specified at t = 0, and continuity is not enforced at the element boundary tn for these profiles. As a result, state profiles remain continuous while control profiles are allowed to be discontinuous at element boundaries. In the direct transcription 5114



Article

∫a

b

⎛1 (x 2 + u 2) dt = h⎜ (x0 2 + u0 2) + (x12 + u12) ⎝2 + (x 2 2 + u 2 2) + ... +

⎞ 1 (xN 2 + uN 2)⎟ ⎠ 2

where the discretization interval h = (b − a)/N is constant. In the second method, the corresponding integral term was replaced by Gaussian quadrature. Figure 2 depicts the dynamic optimization results with both PECE and Gaussian quadrature methods displayed. Moreover, in Table 2 a summary of the objective function values is also shown. These values permit a comparison between both methods of addressing the dynamic optimization of the fractional system. In fact, for α = 1 the Gaussian quadrature approach with only 120 elements attains an optimal objective function value that is very close to the analytical solution of (38). Also, while the plots of the optimal control results shown in Figure 2 provide similar responses, the results in Table 1 present a more accurate comparison between the performance of both the PECE and Gaussian quadrature approaches. Nevertheless, both methods have optimal objective function values for α < 1 that follow the same tendency, so we expect that this comparison remains valid. It is also worth noting that these methods were subjected to multiple grid refinements until we reached a proper step size where no further improvement could be observed in accuracy. 4.2. Isothermal Flash. To describe the dynamic operation of a flash unit we will assume isothermal operating conditions, binary mixture, negligible vapor holdup, and vapor−liquid equilibrium conditions. Accordingly, the liquid phase dynamic mass balance for the light component reads as follows: Fz − Vy − Lx dx = dt M

y= Figure 1. Open-loop system response for different values of the fractional order and keeping the control action at u = −0.3. (a) Predictor-corrector scheme: the discretization step was set to 0.001. (b) Gaussian quadrature: 120 finite elements and 3 internal quadrature points.

s.t.

∫0

1

dx = z + ϕ(x − y) − x dτ

(x 2 + u 2 ) d t

Dα x = −x + u , x(0) = 1

u̇ = x + u ,

x(0) = 1 u(1) = 0

(41)

where τ = (tF/M) is the dimensionless time and ϕ = (V/F) is the vaporization ratio. We also modify this model and rewrite the above equation as a fractional order ODE:

(37)

For α = 1, the solution of this problem can be determined analytically from the boundary value problem: x ̇ = −x + u ,

(40)

where x is the liquid phase mole fraction, F is the flow rate of the feed stream, z is the mole fraction of the feed stream, V is the vapor flow rate, L is the liquid flow rate, y is the vapor phase composition, β is the relative volatility, M is the liquid holdup, and t is the processing time. To avoid scaling and consistency issues,56,57 the above model is cast in dimensionless form:

Dynamic OptimIzation. The dynamic optimization problem to be solved is posed as a transition problem: 1 min Φ(x , u) = x ,u 2

βx 1 + (β − 1)x

(39)

Dα x = z + ϕ(x − y) − x

(42)

and apply m = ⌈α⌉ initial conditions. As discussed in ref 15, this corresponds to the case of a nonideal residence time with various values of α. For α ∈ (1, 2), we are required to specify two initial conditions even when just a single fractional equation is deployed. Because we are mostly interested in computing transitions between steady-states, we normally set the first initial condition (x00) equal to the initial steady-state while the second initial (x10) condition is set to x10 = 0. In fractional order systems, we normally use as initial conditions those conditions specified for an integer order equation. Moreover, additional derivatives that need to be specified are set to zero. A rigorous

(38a) (38b)

and the optimal value of the objective function is Φ* ≈ 0.192 909 3. When solving the above dynamic optimization problem with the PECE algorithm, the integral term in the above dynamic optimization formulation is computed using the trapezoidal rule (the same procedure was deployed in the rest of the examples where PECE was applied): 5115



Article

Figure 2. Optimal system response and control profile for different values of the fractional order for the fundamental linear ordinary differential equation. (a and b) Predictor-corrector scheme: step size = 0.01. (c and d) Gaussian quadrature: 120 finite elements and 3 internal quadrature points.

and different values of the fractional order. The dynamic simulation results obtained using both the PECE and Gaussian quadrature algorithms match very well. It is interesting to note that as the fractional order becomes smaller (α = 0.5) the dynamic response turns out to be slower and more damped, a fact that has been also observed in other processing equipment.46 On the other hand, as the fractional order is raised (α = 1.8), the system response starts to look like a typical second order underdamped system. Eventually, if we keep increasing the fractional order, the system response will blow up like the observed behavior around an unstable steadystate. Because of the nonlinearities embedded in the system, the dynamic response turns out to be asymmetric around the nominal steady-state operating point. This fact can be appreciated from the system response curves for ϕ = 0.2 and ϕ = 0.7. Dynamic Optimization. We now assess the optimal fractional time dynamic behavior of the isothermal flash unit for transitions between the nominal steady-state (point A in Figure 3) and the steady-states denoted by points B and C (see Figure 3). For this case study the dynamic optimization problem was formulated as

Table 1. Linear Fractional ODE System: Objective Function Values for the PECE and Gaussian Quadrature Approaches α

objective function (PECE)

objective function gaussian quadrature

1 0.9375 0.875 0.75 0.5

0.1929175 0.1844517 0.1763457 0.1612493 0.1355693

0.1929062 0.1844351 0.1763210 0.1611913 0.1351318

mathematical proof that justifies this manner of setting initial conditions can be found elsewhere.58 Fractional Open-Loop Dynamic Behavior. We first consider the dynamic behavior of a flash unit for different order of the fractional model. For this we consider the following data for steady-state design: F = 100 g mol/s, z = 0.4, and β = 12.36. Under these conditions, the steady-state values are given as follows: V = 41.856 g mol/s, L = 58.144 g mol/s, x = 0.171, and y = 0.718. Figure 3 depicts the complete steady-state diagram for the flash unit. Hence, we analyze the dynamic response around this operating point. Figure 4 displays the dynamic response for different values of the vaporization ratio 5116



Article

Figure 3. Monotonic steady-state map of the isothermal flash unit. A is the nominal operating design point, whereas B and C are other two operating points with coordinates: (0.2, 0.2955) and (0.7, 0.0868), respectively. tf

min

∫0

s.t.

Dα x = z + ϕ(x − y) − x , x(0) = xAss

x ,ϕ

(x(t ) − x t )2 + (ϕ(t ) − ϕt )2 dt (43)

where xt is the desired end value of the liquid phase composition, ϕt is the reference value of the manipulated variable needed to compute xt, xssA is the steady-state value of the liquid phase mole fraction of the nominal operating point A, and tf denotes the transition time. Figure 5 displays the dynamic optimization results for the A → B and A → C transitions for different fractional order values of the dynamic model and using both the PECE and Gaussian quadrature algorithms for solving the optimal control problem. From Figure 5b, it is interesting to notice that the behavior of the manipulated variable (ϕ) turned out to be a simple step change from the initial to the final value. This probably has to do with the fact that, as seen from Figure 3, the isothermal flash model does not feature highly nonlinear behavior around the addressed operating region and there is only a single steady state. Regarding the quality of the optimal transitions depicted in Figure 5, it seems that both PECE and Gaussian quadrature exhibit similar responses. But, in fact, the responses are not completely equivalent. On one hand, the shape of the control actions is not the same and, as shown in Table 2, the Gaussian quadrature method provides better optimal solutions. Again, this comparison was done with a sufficiently large number of discretization points for both methods. As an additional comparison, this example with α = 1 was solved using direct transcription with three quadrature points and 1000 finite elements for the optimal control problem. For the transitions A → B and A → C, the objective function values were 1.24812 × 10−4 and 2.9679 × 10−5, respectively. 4.3. Bioreactor. In this case study we analyze the model of an isothermal biochemical reactor proposed by Agrawal et al.59 The dimensionless fractional order model reads as follows: Dα1x1 = −x1 + DaM(x 2)x1

(44)

Dα2x 2 = −x 2 + DaΣ(x 2)x1

(45)

Figure 4. Step dynamic response of the isothermal flash unit for different fractional orders and starting from the nominal operating steady-state (point A in Figure 3): α = 0.5 (), α = 1 (- - -), α = 1.5 (· -) and different values of the vaporization ratio (ϕ). (a) Predictorcorrector scheme: step size = 0.3175. (b) Gaussian quadrature: 120 finite elements and 3 interior quadrature points.

where x1 is the normalized cell concentration, x2 is the substrate conversion, Da is the Damköhler number, M(x2) = (1 − x2)ex2/γ is the normalized specific growth rate, and Σ(x2) = M(x2)(1 + β)/(1 + β − x2) is the normalized substrate consumption rate. Setting β = 0.1 and γ = 0.4, the bioreactor depicts the (steadystate) multiplicity pattern shown in Figure 6. Also, α1 and α2 are the fractional orders of the first and second differential equations, respectively. As discussed in ref 15, the fractional derivatives in this model correspond to nonlocal diffusion and nonideal residence times, which are typical of cellular reactors and suspended particle systems. Note that the effect of these phenomena will generally be different for each state, the cell concentration and substrate conversion. Fractional Open-Loop Dynamic Behavior. The dynamic simulation results for the bioreactor case study using both the PECE and Gaussian quadrature approaches are compared in Figure 7. Here we notice that as the Damköhler number decreases the reactor dynamic response becomes oscillatory, indicating the proximity of the unstable region (as seen from Figure 6) and a Hopf bifurcation point. However, the dynamic response of the fractional order system remains nonoscillatory. 5117



Article

Figure 5. Dynamic optimization results for the isothermal flash unit for transitions from the nominal operating steady-state (A) to the steady states B and C for different values of the fractional order. (a and b) Predictor-corrector scheme: step size = 0.1508. (c and d) Gaussian quadrature: 120 finite elements and 3 interior quadrature points.

Table 2. Flash System: Objective Function Values for the PECE and Gaussian Quadrature Approaches transition

α

A→B

1.5 1 0.5 1.5 1 0.5

A→C

objective function PECE 1.5004 1.1364 2.1739 9.3542 8.0921 8.9106

× × × × × ×

10−2 10−2 10−2 10−3 10−3 10−3

objective function Gaussian quadrature 1.8718 1.2414 2.8409 5.4324 2.9316 3.5540

× × × × × ×

10−4 10−4 10−4 10−5 10−5 10−5

A way to explain this behavior for the integer system lies in the fact that as the Damköhler number of the system is decreased, the magnitude of the imaginary part of the eigenvalues of the equivalent linear system increases, leading to the observed oscillatory behavior. It is clear that for the same change in the Damköhler number, the dynamic behavior of the corresponding fractional order system turns out to be faster and nonoscillatory. Therefore, we can conclude that the magnitude of the absolute value of the eigenvalues of the fractional order system increases, whereas the magnitude of their imaginary parts does not increase as much.

Figure 6. Multiple steady-state map of the bioreactor system. The continuous line stands for open-loop stable steady states while the dashed line represents open-loop unstable steady-states.

5118



Article

Figure 7. (a and c) Step open-loop system response for the fractional order bioreactor system (α1 = α2 = 0.8; solid line) and first-order models (dashed line) for differentes values of the Damkholer number. (b and d) Step open-loop system response for different fractional orders and Da = 0.7. (a and b) Predictor-corrector scheme: step size = 0.001, (c) Gaussian quadrature: 120 finite elements and 3 internal quadrature points, (d) Gaussian quadrature: 1200 finite elements and 3 internal quadrature points.

In summary, fractional order models with α < 1 seem to give rise to smoother time responses. Commonly, from a control point of view, we would like to avoid oscillatory behavior, so dynamic fractional order models seem to offer this advantage. However, when it comes to analyzing the bioreactor dynamic response for fractional orders greater than the nominal one the situation turns out to be different as shown in Figures 7b and d. In these figures, the advantage of using a fractional order system smaller than the nominal one are more evident: for a slight increase of the nominal order (α1 = α2 = 1.2) the dynamic system response depicts sustained oscillatory behavior which is a clear indication of a Hopf bifurcation point. In fact, with an additional small fractional order increase the system would turn out to be unstable. Finally, as seen from Figure 7 the system response using the PECE and Gaussian quadrature algorithms matches well. Dynamic Optimization. We now conduct dynamic optimization calculations for α < 1, as these lead to nonoscillatory behavior. Here we examine the fractional transition dynamic optimal behavior from the nominal operating point A (see Figure 6) to the steady-states denoted as B, C, and D. As noticed

from Figure 6, all the final steady-states are open-loop unstable. The dynamic optimization problem is given by tf

x1, x 2 , u

∫0

s.t.

Dα1x1 = −x1 + DaM(x 2)x1 , x1(0) = x1ss

min

((x1(t ) − x1t )2 + (x 2(t ) − x 2t )2

+ (Da(t ) − Dat )2 dt Dα2x 2 = −x 2 + Da Σ(x 2)x1 , x 2(0) = x 2ss

(46)

where the superscript t stands for the processing conditions to be reached after system transition, and the superscript ss represents the initial steady-state conditions. The manipulated variable is the Damköhler number Da. In all the addressed cases, the fractional order of the underlying system was changed as follows: (a) The fractional order of the first equation (α1) was set to one, whereas the fractional order of second equation (α2) took the following values: 1, 0.8, and 0.6. (b) This time α2 = 1 and α1 took the values 0.8 and 0.6. (c) Finally, the order of both equations was simultaneously changed and set to the same value α1 = α2 = 0.8 and 0.6. The aim of running these calculations was to get a broad picture of the fractional order 5119



Article

Figure 8. PECE: dynamic optimization results of the bioreactor system for the A → B transition. (a and b) The order of the first ODE was set to 1, and the order of the second ode was changed. (c and d) The order of the second ODE was set to 1, and the order of the first ode was changed. (e and f) The fractional order of the two ODEs was simultaneously changed according to the displayed values.

For the A → B transition from Figures 8 and 9, we see that the dynamic response of the bioreactor is satisfactory in almost all the addressed cases. This turns out to be the case especially for Figure 8a. Regarding Figure 8c and e, there is a small change in the state variables at the end of the transition period. It turns

bioreactor system under optimal control conditions. We also remark that for the PECE algorithm the step size was kept at 0.0224, and for the Gaussian quadrature method, 240 elements finite with 3 internal quadrature points were deployed. 5120



Article

Figure 9. Gaussian quadrature: dynamic optimization results of the bioreactor system for the A → B transition. (a and b) The order of the first ODE was set to 1, and the order of the second ode was changed. (c and d) The order of the second ODE was set to 1, and the order of the first ode was changed. (e and f) The fractional order of the two ODEs was simultaneously changed according to the displayed values.

out that both PECE and Gaussian quadrature approaches give essentially the same results with very small differences at the end of the transition periods. In the case of the A → C transition, the results are shown in Figures 10 and 11. These results indicate that for fractional orders α2 = 0.8 and 0.6 the

performance of optimal control system outperforms that obtained from using integer order dynamic models, in the sense that smaller control efforts seem to be necessary to reach the target steady-state. This feature can be highlighted as an 5121



Article

Figure 10. PECE: dynamic optimization results of the bioreactor system for the A → C transition. (a and b) The order of the first ODE was set to 1, and the order of the second ode was changed. (c and d) The order of the second ODE was set to 1, and the order of the first ode was changed. (e and f) The fractional order of the two ODEs was simultaneously changed according to the displayed values.

advantage of the fractional system in comparison to the integer order system. For transition A → D both PECE and Gaussian quadrature approaches render similar results as seen from Figures 12 and 13. However, the quality of the response of the fractional order

dynamic system gets worse as the order of the plant decreases as shown in part (e) of both figures. This is especially true for the situation depicted in Figure 12f. This situation can be partially fixed by adding weighting functions to the corresponding objective functions. For this transition, the best dynamic 5122



Article

Figure 11. Gaussian quadrature: dynamic optimization results of the bioreactor system for the A → C transition. (a and b) The order of the first ODE was set to 1, and the order of the second ode was changed. (c and d) The order of the second ODE was set to 1, and the order of the first ode was changed. (e and f) The fractional order of the two ODEs was simultaneously changed according to the displayed values.

optimal results were obtained by fixing α = 1 for the second ODE of the underlying model and setting the order of the first ODE to 0.8 as shown in parts c and d of Figures 12 and 13. In summary, we can conclude that, in general terms, for the

fractional order bioreactor system the amount of control action turns out to be smaller than when considering the same plant with integer order. As the size of the intended transition increases, better control actions were attained using fractional 5123



Article

Figure 12. PECE: dynamic optimization results of the bioreactor system for the A → D transition. (a and b) The order of the first ODE was set to 1, and the order of the second ODE was changed. (c and d) The order of the second ode was set to 1, and the order of the first ode was changed. (e and f) The fractional order of the two ODEs was simultaneously changed according to the displayed values.

orders closer to the nominal one. However, for small size transitions, smaller fractional orders give rise to improved control actions. Although the optimal dynamic responses displayed in Figures 10−13 look similar, from Table 3 we see

that broadly the Gaussian quadrature method provides better optimal solutions than the PECE method. As an additional comparison, we solved this example with α1 = α2 = 1 using direct transcription with three quadrature points 5124



Article

Figure 13. Gaussian quadrature: dynamic optimization results of the bioreactor system for the A → D transition. (a and b) The order of the first ODE was set to 1, and the order of the second ODE was changed. (c and d) The order of the second ode was set to 1, and the order of the first ode was changed. (e and f) The fractional order of the two ODEs was simultaneously changed according to the displayed values.

and 1000 finite elements for the optimal control problem. For the transitions A → B, A → C, and A → D, the objective function values were 0.0975097, 0.0273253, and 0.486435, respectively.

5. CONCLUSIONS In this work we address both the dynamic simulation and optimization of mathematical models that can be represented in terms of fractional differential equations (FDEs). Although the 5125



Article

Table 3. Bioreactor System: Objective Function Values for the PECE and Gaussian Quadrature Approaches PECE

Gaussian quadrature

α1

α2

A→B

A→C

A→D

A→B

A→C

A→D

1 1 1 0.8 0.6 0.8 0.6

1 0.8 0.6 1 1 0.8 0.6

0.076789 0.061684 0.051033 0.074571 0.075365 0.061827 0.055184

0.204235 0.195573 0.207907 0.201556 0.200103 0.193329 0.207099

0.505019 0.559829 0.657928 0.513845 0.524236 0.579308 0.721828

0.09750 0.10067 0.10547 0.09836 0.09914 0.10136 0.10588

0.02732 0.02844 0.03195 0.02692 0.02623 0.02800 0.03048

0.04864 0.05445 0.06553 0.04905 0.04901 0.05577 0.06926

Notes

concept of fractional derivatives goes back almost two centures, very few studies deal with the simulation and optimal control of FDE systems. This is partly due to the lack of general numerical algorithms for these dynamic systems. Existing algorithms are often tailored to deal with specific applications and are generally intended for small scale systems. For the applications in this study, we have deployed a predictor-corrector method proposed in the literature22 and we have developed a novel Gaussian quadrature method to address such systems. While both methods can be used for dynamic simulation and optimization purposes, our results indicate the superiority of the Gaussian quadrature approach. This is expected as the predictor-corrector method is an explicit second order method, while the higher order Gaussian quadrature approach features better approximation properties. This behavior is also seen for FDEs in the addressed case studies; as the effects of nonlinear dynamics increase the performance of the PECE method gets worse, as accuracy is harder to maintain. However, the corresponding performance of the proposed Gaussian quadrature method maintains its approximation level even in the face of complex dynamic behavior. These results encourage us to apply this method to more complex fractional order systems. Regarding the precision of the Gaussian quadrature approach when compared against the PECE method we would like to highlight that the precision of each method also depends upon the size of the discretization mesh. We refined the results from each method using the proper number of discretization points. Starting from an arbitrary number of discretization points, we just kept increasing/decreasing the number of discretization points until the response of the system remains constant. This is the reason why sometimes the number of discretization points used in each method is not exactly the same. Detailed stability and accuracy properties of the FDE-based Gaussian quadrature method still remain to be analyzed. Our results indicate that in some cases smaller and more damped control actions could be demanded if the model is formulated in terms of a fractional system. In future work we will explore the use of larger and more complex fractional models, especially in modeling new energy sources such as batteries and fuel cells. Moreover, we think that fractional models extend very well in modeling convection and dispersion in the mixing of chemical reactors or in trayed distillation columns. In addition, areas such as system biology and nanotechnology with complex geometries, which could be modeled with fractal dimensions, would benefit greatly from the deployment of fractional order systems.

■

The authors declare no competing financial interest.

■

DEDICATION We are honored to dedicate this study to the memory of Profs. David Himmelblau and Gary Powers. Both were pioneers in Process Systems Engineering, especially in pursuing challenging research topics outside of the mainstream. It is in this spirit that we prepared this work on fractional differential equations.

■

REFERENCES

(1) Oldham, K. B.; Spanier, J. The Fractional Calculus; Academic Press, 1974. (2) Miller, K. S.; Ross, B. An Introduction to the Fractional Calculus and Fractional Differential Equations; John Wiley, 1993. (3) Podlubny, I. Fractional Differential Equations; Academic Press, 1999. (4) Hilfer, R. Applications of Fractional Calculus in Physics; World Scientific, 2000. (5) Bonilla, B.; Rivero, M.; Rodriguez-Germa, L.; Trujillo, J. Fractional differential equations as alternative models to nonlinear differential equations. Appl. Math. Comput. 2007, 187, 79−88. (6) Koeller, R. C. Applications of fractional calculus to the theory of viscoelasticity. J. Appl. Mech. 1984, 51 (2), 229−307. (7) Momani, S. An algorithm for solving the fractional convectiondiffusion equation with nonlinear source term. Communications in Nonlinear Science and Numerical Simulation 2007, 12, 1283−1290. (8) Ochoa-Tapia, J. A.; Valdes-Parada, F. J.; Alvarez-Ramirez, J. A fractional-order darcy’s law. Phys. A 2007, 374 (1), 1−14. (9) Axtell, M.; Bise, M. E. Fractional calculus applications in control systems. In Proceedings of the National Aerospace and Electronics Conference, Dayton, Ohio, May 21−25; IEEE, 1990; pp 563−566. (10) Laskin, N. Fractional market dynamics. Phys. A 2000, 287, 482− 492. (11) Sokolov, I. M.; Klafter, J.; Blumen, A. Fractional kinetics. Phys. Today 2002, No. November, 48−54. (12) Dokoumetzidis, A.; Macheras, P. Fractional kinetics in drug absorption and disposition processes. J. Pharmacokinet. Pharmacodyn. 2009, 36, 165−178. (13) Ionescu, C. Fractal structure and storage dynamics of glycogen. Proceedings of the 9th International Symposium on Dynamics and Control of Process Systems (DYCOPS 2010), Neuven, Belgium, July 5−7; IFAC, 2010; pp 330−334. (14) Machado, J. A.; Costa, A. C.; Qhelhas, M. Fractional dynamics in DNA. Commun. Nonlin. Sci. Numer. Simulat. 2011, 16, 2963−2969. (15) Schumer, R.; Meerschaert, M.; Baeumer, B. Fractional advectiondispersion equations for modeling transport at the Earth surface. J. Geophys. Res 2009, 114, F00A07 DOI: 10.1029/2008JF001246. (16) Atici, F. M.; Sengul, S. Modeling with fractional difference equations. J. Math. Anal. Appl. 2010, 369 (1), 1−9. (17) Song, L.; Xu, S.; Yang, J. Dynamic models of happiness with fractional order. Commun. Nonlin. Sci. Numer. Simulat. 2010, 15, 616− 628.

AUTHOR INFORMATION

Corresponding Author

*On leave from Universidad Iberoamericana, México. E-mail: antonio.fl[email protected]. 5126



Article

(43) Samko, S. A.; Kilbas, A. A.; Marichev, O. I. Fractional Integrals and Derivatives: Theory and Applications; Gordon and Breach, 1993. (44) Predictor-corrector PECE method for fractional differential equations. http://www.mathworks.com/matlabcentral/fileexchange/ 32918 (accesed on Feb 7, 2011). (45) Biegler, L. T. Nonlinear Programming: Concepts, Algorithms and Applications to Chemical Engineering; SIAM, 2010. (46) Ahmad, W.; Abdel-Jabbar, N. Modeling and simulation of a fractional order bioreactor system. Workshop on Fractional Differentiation and its Applications, Portugal, July 19−21; IFAC, 2006. (47) Sabatier, J.; Cugnet, M.; Laruelle, S.; Grugeon, S.; Sahut, B.; Oustaloup, A.; Tarascon, J. M. A fractional order model for lead-acid battery crankability estimation. Commun. Nonlin. Sci. Numer. Simulat. 2010, 15 (5), 1308−1317. (48) Cao, H.; Deng, Z.; Li, X.; Yang, J.; Qin, Y. Dynamic modeling of electrical characteristics of solid oxide fuel cells using fractional derivatives. Int. J. Hydrogen Energy 2010, 35 (4), 1749−1758. (49) Bertrand, N.; Sabatier, J.; Briat, O.; Vinassa, J. Fractional nonlinear modelling of ultracapacitors. Commun. Nonlin. Sci. Numer. Simulat. 2010, 15, 1327−1337. (50) Denga, Z.; Caoa, H.; Lia, X.; Jianga, J.; Yang, J.; Qin, Y. Generalized predictive control for fractional order dynamic model of solid oxide fuel cell output power. J. Power Sources 2010, 195, 8097− 8103. (51) Dulf, E.; Both, R.; Dumitrache, D. C. Fractional order models for a cryogenic separation column. Proceedings of the International Conference on Automation Quality and Testing Robotics, Romania, May 28−30; IEEE, 2010; pp 1−6. (52) Fourer, R.; Gay, D. M.; Kernighan, B. W. AMPL: A Modeling Language for Mathematical Programming, 2nd ed.; Duxbury Press: Pacific Grove, CA, 2003. (53) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R. A. GAMS: A User’s Guide; GAMS Development Corporation, 1998; http://www. gams.com. (54) Waechter, A.; Biegler, L. T. On the implementation of an interior point filter line search algorithm for large-scale nonlinear programming. Math. Progr. 2006, 106 (1), 25−57. (55) de Hoog, F. R.; Knight, J. H.; Stokes, A. N. An improved method for numerical inversion of laplace transforms. SIAM J. Sci. Stat. Comp. 1982, 3 (3), 357−366. (56) Popovic, J. K.; Atanackovic, M. T.; Pilipovic, A. S.; Rapaic, M. R.; Pilipovic, S.; Atanackovic, T. M. A new approach to the compartmental analysis in pharmacokinetics: fractional time evolution of diclofenac. J. Pharmacokinet. Pharmacodyn. 2010, 37 (2), 119−134. (57) Dokoumetzidis, A.; Magin, R.; Macheras, P. A commentary on fractionalization of multi-compartmental models. J. Pharmacokinet. Pharmacodyn. 2010, 37 (2), 203−207. (58) Diethelm, K.; Ford, N. J. Numerical solution of the bagley torvik equation. BIT Numer. Math. 2002, 42 (3), 490−507. (59) Agrawal, P.; Lee, C.; Lim, H. C.; Ramkrishna, D. Theoretical Investigations of Dynamic Behav- ior of Isothermal Continuous Stirred Tank Biological Reactors. Chem. Eng. Sci. 1982, 37 (3), 453−462. (60) (a) Samko, S. G.; Kilbas, A. A.; Marichev, O. I. Fractional Integrals and Derivatives Theory and Applications; Gordon and Breach: New York, 1993. (b) Oldham, K. B.; Spanier, J. The fractional Calculus; Academic Press: New York, 1974.

(18) Du, Q.; Gunzburger, M.; Lehoucq, R. B.; Zhou, K. Analysis and approximation of nonlocal diffusion problems with volume constraints. SIAM Rev. 2013, 54 (4), 667−696. (19) Diethelm, K.; Ford, N. J.; Freed, A. D.; Luchko, Y. Algorithms for the fractional calculus: A selection of numerical mathods. Comput. Methods Appl. Mech. Eng. 2005, 194 (6−8), 743−773. (20) Griffiths, D. F.; Higham, D. J. Numerical Methods for Ordinary Differential Equations; Springer-Verlag, 2010. (21) Ascher, U. M.; Petzold, L. R. Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations; SIAM, 1998. (22) Diethelm, K.; Ford, N. J.; Freed, A. D. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. 2002, 29 (1−4), 3−22. (23) Rawashdeh, E. A. Numerical solution of fractional integrodifferential equations by collocation method. Appl. Math. Comput. 2006, 176 (1), 1−6. (24) Diethelm, K. The Analysis of Fractional Differential Equations; Springer-Verlag, 2010. (25) McAuley, K. B.; MacGregor, J. F. Optimal grade transitions in a gas phase polyethylene reactor. AIChE J. 1992, 38 (10), 1564−1576. (26) Flores-Tlacuahuac, A.; Biegler, L. T.; Saldivar-Guerra, E. Optimization of HIPS Open-Loop Unstable Polymerization Reactors. Ind. Eng. Chem. Res. 2005, 44, 2659−2674. (27) Bansal, V.; Perkins, J. D.; Pistikopoulos, E. A Case Study in Simultaneous Design and Control Using Rigurous, Mixed-Integer Dynamic Optimization Models. Ind. Eng. Chem. Res. 2002, 41, 760− 778. (28) Sakizlis, V.; Perkins, J. D.; Pistikopoulos, E. Recent Advandces in Optimization-Based Simultaneous Process and Control Design. Comput. Chem. Eng. 2004, 28, 2069−2086. (29) Zavala, V. M.; Biegler, L. T. Optimization-based strategies for the operation of low-density polyethylene tubular reactors: Nonlinear model predictive control. Comput. Chem. Eng. 2009, 33 (10), 1735− 1746. (30) Agrawal, O. M. Formulation of Euler-Lagrange equations for fractional variational problems. J. Math. Anal. Appl 2002, 272 (1), 368−379. (31) Agrawal, O. M. A general formulation and solution scheme for fractional optimal control problems. Nonlin. Dynam. 2004, 38 (1−4), 323−337. (32) Agrawal, O. M. A quadratic numerical scheme for fractional optimal control problem Journal of Dynamic Systems. Measurement Control 2008, 130 (1), 1−6. (33) Agrawal, O. M. A formulation and numerical scheme for fractional optimal control problems. J. Vibration Control 2008, 149 (9− 10), 1291−1299. (34) Baleanu, D.; Defterli, O.; Agrawal, O. M. A central diefference numerical scheme fo fractional optimal control problems. J. Vibration Control 2009, 15 (4), 583−597. (35) Tricaud, C.; Chen, Y. Time-optimal control of systems with fractional dynamics. Int. J. Differential Eq. 2010, No. 461048, DOI: 10.1155/2010/461048. (36) Bryson. A. E.; Ho, Y. C. Applied Optimal Control; Taylor and Francis, 1981. (37) Biegler, L. T. Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes; SIAM, 2010. (38) Diethelm, K. Multi-term fractional differential equations, multiorder fractional differential systems and their numerical solution. J. Eur. Syst. Automatiss 2008, 28 (6−8), 665−676. (39) Verotta, D. Fractional dynamics pharmacokinetics pharmacodynamic models. J. Pharmacokinet. Pharmacodyn. 2010, 37, 257−276. (40) Folland, G. B. Advanced Calculus; Prentice Hall, Upper Saddle River, NJ, 2002. (41) Caputo, M. Linear models of dissipation whose Q is almost frequency independent. Geophys. J. R. Astronom. Soc. 1967, 13, 529− 539. (42) Diethelm, K.; Ford, N. Analysis of fractional diferential equations. J. Math. Anal. Appl. 2002, 265 (2), 229−248. 5127


Optimization of Fractional Order Dynamic Chemical Processing Systems

Recommend Documents