Differential-geometric model-based control (DGMBC): A software package for controller design

Available online at www.sciencedirect.com

Computers and Chemical Engineering 32 (2008) 1569–1588

Differential-geometric model-based control (DGMBC): A software package for controller design Chanin Panjapornpon a , Masoud Soroush a,∗ , Warren D. Seider b b

a Department of Chemical and Biological Engineering, Drexel University, Philadelphia, PA 19104, United States Department of Chemical and Biomolecular Engineering, University of Pennsylvania, Philadelphia, PA 19104-6393, United States

Received 22 August 2005; received in revised form 23 July 2007; accepted 26 July 2007 Available online 2 August 2007

Abstract This paper presents a new software package, called differential-geometric model-based control (DGMBC), which carries out symbolic manipulations to automatically generate differential-geometric, model-based controllers and subsequently tests the designed controller. This prototype software was developed to simplify the industrial implementation and testing of differential-geometric, model-based controllers on lumpedparameter processes. DGMBC has a user-friendly interface that allows a user to enter process model equations and parameters easily. The user interface was developed using Visual Basic and linked to MATHEMATICA using MathLink. The user enters the process model (set of ordinary differential equations), and the software generates an analytical model-based controller, if an analytical solution exists. The resulting analytical model-based controller (set of ordinary differential and algebraic equations) can be in FORTRAN, C, or MATLAB format. DGMBC can also simulate the closed-loop process responses. The application and implementation of DGMBC 1.0 are shown using three chemical and biochemical process examples with varying levels of complexity. An analytical model-based controller is designed for each of the processes, and simulation results showing the closed-loop process responses are presented. © 2007 Elsevier Ltd. All rights reserved. Keywords: Controller design software; Model-based control; Input–output linearization; Feedback linearization; Numerical simulation; Computer-aided design software

1. Introduction Model-predictive control and differential-geometric control have been the most widely used model-based control methods. In model-predictive control, the controller action is the solution of an optimization problem that is solved numerically on-line. In contrast, differential-geometric control is a direct synthesis approach in which the controller is derived by requesting a desired closed-loop response in the absence of input constraints. Both control methods lead to controllers that include an inverse of a model of the process under consideration. In model-predictive control, the model inversion is numerical, while in differential-geometric control, it is analytical. The analytical inversion in differential-geometric control requires analytical partial derivatives and symbolic manipulations, which may become cumbersome as the relative order and/or the level of complexity of the model increases. Indeed, the burden of taking analytical partial derivatives and performing symbolic manipulations has been a major factor preventing this theoretically sound, efficient control method from being implemented widely in the process industries. The increasing role of computers in solving complex scientific and engineering problems has led to the development of software packages that facilitate the design of model-based controllers (Barker, Grant, & Zhuang, 1996; Barker & Zhuang, 1997; Bemporad & Morari, 2001; Blankenship, Kwatny, LaVigna, & Polyakov, 1997; Campbell, Delebecque, & von Wissel, 1996; Datta & Sarkissian, 1999; de Jager, 1995; Kitamoto, Saeki, & Ando, 1992; Ohtani, Fukuzawa, & Masubuchi, 1994; Sack & Singh, 2000). Software such as MATHEMATICA (Mathematica, 1999), MAPLE (Maple, 2000), and the Symbolic Toolbox in MATLAB (Matlab, 2002) have ∗

Corresponding author. Tel.: +1 215 895 1710; fax: +1 215 895 5837. E-mail address: [email protected] (M. Soroush).

0098-1354/$ – see front matter © 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2007.07.006

1570

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

Nomenclature A, B, C chemical species cp heat capacity of feed and product (kJ kg−1 K−1 ) cpj heat capacity of jacket fluid (kJ kg−1 K−1 ) CA0 inlet concentration of the reactant A (mol l−1 ) CA , CA1 , CA2 outlet concentration of the reactant A (mol l−1 ) CB , CB1 , CB2 concentration of B (mol l−1 ) D differential operator (D = d/dt) Ea , E1 , E2 activation energy (kJ kmol−1 ) F reactor flow rate (m3 h−1 or l h−1 ) jacket coolant flow rate (l h−1 ) Fj −H1 , −H2 heat of reaction (kJ mol−1 ) open loop Jacobian Jol Jcl closed-loop Jacobian Jzd Jacobian of the zero dynamics k0 , k1 , k2 reaction rate constant (h−1 ), (l mol−1 h−1 ), (h−1 ) L matrix of observer gain m number of manipulated inputs n number of controlled outputs p vector of the orders of the requested responses q cooling/heating rate (kJ s−1 ) ri relative order of output variable yi R universal gas constant (kJ kmol−1 K−1 ) S heat transfer surface area (m2 ) t time (s) T reactor outlet temperature (K) Ti reactor inlet temperature (K) Tj jacket temperature (K) Tji jacket inlet temperature (K) u vector of manipulated inputs U overall heat-transfer coefficient (kJ m−2 K−1 s−1 ) V, V1 , V2 reactor volume (m3 ) jacket volume (l) Vj x vector of state variables y vector of controlled outputs ysp vector of set-points Greek letters β adjustable parameter of stability constraint ε1 , . . ., εn adjustable parameters of controller γ1 bioreactor model parameter η state for integral action λ1 ,λ2 ,. . . eigenvalues of Jacobian ρ density of mixture in reactor (kg l−1 ) ρj density of jacket fluid (kg l−1 ) Subscripts A, B chemical species ss steady-state sp set-point 0 initial value, process

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

1571

been used to perform symbolic manipulations. Available software packages can be divided into two categories. In one category, the software is designed as a toolbox for an existing symbolic manipulation package (Blankenship et al., 1997; Datta & Sarkissian, 1999; Kitamoto et al., 1992; Ohtani et al., 1994). For example, Blankenship et al. (1997) developed a package for modeling and designing nonlinear adaptive controllers using MATHEMATICA in which the controller equations are generated in the format suitable for simulation in SIMULINK. However, a user has to program manually sequences of input commands with specific patterns to complete the controller design without an interactive visual interface, and the package cannot be used directly for the controller performance testing and implementation. The software developed by Campbell et al. (1996) also falls in this category. It uses MAPLE to perform the symbolic computations needed for feedback linearization in mechanical multi-body systems, and it then generates a controller code that is suitable for Scilab (www.scilab.com) to carry out closed-loop numerical simulations. Although DGMBC uses an interface between the symbolic manipulation and numerical calculation as in the software packages of Blankenship et al. (1997) and Campbell et al. (1996), the controller design methods that it automates and the interface that it uses are different. In DGMBC, the controller design and testing can be completed within the developed visual interface; DGMBC does not require any additional command codes or numerical simulating software to perform the controller testing. While DGMBC designs model-based controllers based on input–output linearization and shortest-prediction-horizon model-predictive control, the software packages of Blankenship et al. (1997) and Campbell et al. (1996) design feedback linearizing and nonlinear adaptive controllers, respectively. In the second category, multiple calculating engines are used to perform controller design (Campbell et al., 1996; Kitamoto et al., 1992). For example, Kitamoto et al. (1992) created H∞ controller design software in which the design algorithm is programmed in MATHEMATICA and the Front-End interface is written in MATLAB. A user enters the control system in the form of a block diagram. The software is unable to calculate and show the closed-loop performance of the designed controller. Although symbolic computation software packages for model-based controller design have been available for many years, their applications have been limited because of the inadequacy of the available computer-aided design tools and user-interface platforms. For example, the packages do not perform the complete controller design when the process model is in the form of ordinary differential equations. The existing packages require many computational engines to complete the controller design and testing. Furthermore, they are not user-friendly for design integration and closed-loop simulation. Thus, a software package devoid of these limitations should be useful in the practice of process control. Motivated by the deficiencies of the existing software packages for the design of analytical, model-based controllers, in this work, a new software package, called differential-geometric model-based control (DGMBC), is presented. The version presented here is 1.0. DGMBC has a user-friendly interface and fully automates the design of differential-geometric controllers developed for general nonlinear processes. DGMBC generates controller equations in C, FORTRAN, and MATLAB formats and carries out closed-loop simulations. Furthermore, it is able to perform the controller design and testing with the use of only one computational engine. The MATHEMATICA program was selected as the symbolic computational engine, because it has built-in functions that perform sufficiently accurate computations, allowing for the creation of specific programming packages. In DGMBC, a user can develop his/her own functions that support their individual calculation needs on the basis of the available functions in MATHEMATICA. Furthermore, DGMBC does not require the declaration of symbolic variables in an expression, and it supports the communication between an external program (Front-End interface) and its kernel. The Visual Basic language was used for the creation of a stand-alone application that communicates with the MATHEMATICA kernel to perform symbolic manipulations using the MathLink program. The differential-geometric control methods presented in (Daoutidis, Soroush, & Kravaris, 1990; Kanter, Soroush, & Seider, 2002; Kravaris & Soroush, 1990; Panjapornpon & Soroush, 2006; Panjapornpon, Soroush, & Seider, 2006) were programmed as MATHEMATICA packages, to design controllers for general nonlinear processes whether stable, unstable, minimum-phase, and/or non-minimum-phase. DGMBC automates the design of differential-geometric model-based controllers for general nonlinear processes, avoiding laborious analytical calculations and manipulations when the process model is complex and/or has high relative order(s). DGMBC 1.0 has been tested successfully on process models with up to 10 state variables. DGMBC has applications in research, development, training, and education. The software permits the routine design and use of differential-geometric, model-based controllers in the process industries, academia, and research institutions. It enables researchers, practitioners, and instructors to design, implement and study these controllers with ease for complex processes. The most important application of the software is in development; it simplifies greatly the design, implementation and testing of differential-geometric, model-based controllers on industrial processes. The software can also be used in research studies on nonlinear control and process design analyses, especially when the process is non-minimum-phase; an example is the study presented in (Meel, Seider, & Soroush, 2006). In addition to the development and research applications, the software can be used as a teaching aid in personnel training and in nonlinear control and/or analysis courses. The paper is organized as follows. The scope of the study and some mathematical preliminaries on differential-geometric, modelbased controller design which is automated by DGMBC, are given in Section 2. Section 3 describes several aspects of DGMBC 1.0. Finally, the application and implementation of DGMBC 1.0 are illustrated in Section 4. Concluding remarks are given in Section 5.

1572

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

2. Scope and mathematical preliminaries DGMBC 1.0 synthesizes the differential-geometric control methods presented in (Daoutidis et al., 1990; Kanter et al., 2002; Kravaris & Soroush, 1990; Panjapornpon & Soroush, 2006; Panjapornpon et al., 2006). The synthesized controller always includes a state feedback and a dynamic system to add integral action. It may include a state observer to reconstruct unmeasured state variables. In this section, the state feedback design is reviewed briefly. For the observer design, the reader can refer to (Soroush, 1998). DGMBC 1.0 is applicable to the general class of multivariable processes having a mathematical model in the form: dx dt y

=

f (x, d, u),

=

h(x)

x(0) = x0

(1)

where x = [x1 · · ·xn ]T ∈ Rn is the vector of state variables, u = [u1 · · ·um ]T ∈ Rm is the vector of manipulated inputs, d = [d1 · · ·dq ]T ∈ Rq is the vector of measured disturbances, y = [y1 · · ·ym ]T ∈ Rm is the vector of controlled outputs, and f(x,d,u) = [f1 (x,d,u)· · ·fn (x,d,u)]T and h(x) = [h1 (x)· · ·hm (x)]T are smooth functions. For the class of processes described by (1), the following assumptions are made: • Relative orders (degrees) of the process outputs, y1 , . . ., ym , respectively denoted by r1 , . . ., rm , are locally finite and constant. The relative order, ri , is the smallest integer for which dri yi /dt ri depends explicitly on u. • The characteristic (decoupling) matrix of the process ⎤  ∂ dr1 y1 ⎢ ∂u dt r1 ⎥ ⎢ ⎥ ⎢ ⎥ .. ⎢ ⎥ . ⎢ ⎥ ⎢ ⎥  ⎣ ∂ drm ym ⎦ dt rm ∂u ⎡

is locally nonsingular. This assumption is needed only for the two control methods described in Subsections 2.1 and 2.3.2. • The process is locally controllable and observable. • The Jacobian of the process, ∂f(x,d,u)/∂x, and [∂h(x)/∂x][∂f(x,d,u)/∂x]−1 [∂f(x,d,u)/∂u] evaluated at the desired steady-state are nonsingular. The following notations are used for each hi , i = 1, . . ., m: h0i (x)  hi (x)

−1 ˙ ¨ ˙ d, ¨ . . .)  ∂hi (x, d, d, d, . . .) f (x, d, u), = 1, . . . , ri − 1 h i (x, d, d, ∂x ri −1 (x, d, d, ˙ d, ¨ . . .) ∂h i r ˙ d, ¨ · · ·)  hii (x, u, d, d, f (x, d, u) ∂x 2.1. Minimum-phase processes For this sub-class of processes, one can implement the input–output linearization method described in (Kravaris & Soroush, 1990) for processes without measured disturbances and in (Daoutidis et al., 1990) for processes with measured disturbances. For a process in the form of (1), closed-loop output responses of the following form are requested: (ε1 D + 1)r1 y1 = ysp1 .. .

(2)

(εm D + 1) ym = yspm rm

where D = d/dt, and ε1 , . . ., εm are the time constants of the process outputs, y1 , . . ., ym , in closed-loop, respectively. Using the model of (1), the time derivatives of the outputs are calculated analytically, which requires taking partial derivatives with respect to x and performing matrix operations. After substituting for the time derivatives in (2), a feedforward/state feedback of the following form is calculated: ˙ d, ¨ . . .) u = ψr (x, ysp , d, d,

(3)

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

1573

2.2. Stable non-minimum-phase processes Kanter et al. (2002) extended the method in (Kravaris & Soroush, 1990) to stable, non-minimum-phase, nonlinear processes, leading to an approximate input–output linearizing control method. For the processes in the form of (1), linear, closed-loop, processoutput responses of the following form are requested: (ε1 D + 1)p1 y1 = ysp1 .. . pm (εm D + 1) ym = yspm

(4)

where p1 ≥ r1 , . . ., pm ≥ rm . Using the model of (1), the time derivatives of the outputs are calculated analytically, which again requires taking partial derivatives with respect to x and performing matrix operations. Note that the time derivatives of an output yi that are of order higher than ri depend also on the time derivative(s) of u. After substituting for the time derivatives of the controlled outputs in (4), one obtains a feedforward/dynamic state feedback in the compact form: ˙ d, ¨ . . .) = ysp p (x, u, u(1) , u(2) , . . . , d, d,

(5)

After setting all time derivatives of u in (5) to zero, if (5) is solvable for u, a feedforward/static state feedback in the general form: ˙ d, ¨ . . .) u = ψp (x, ysp , d, d,

(6)

is calculated. The tunable parameters ε1 , . . ., εm are chosen such that for every l = 1, . . ., m, all of the eigenvalues of [I + εl Jol (xss ,dss ,uss )] lie inside the unit circle, where Jol (xss ,dss ,uss ) is the Jacobian of (1) evaluated at the desired steady-state. 2.3. General processes The application of the method described in (Kanter et al., 2002) is limited to stable processes. To design controllers that are also applicable to unstable, non-minimum-phase processes, Panjapornpon et al. (2006) and Panjapornpon and Soroush (2006) proposed the two controller design methods described below. 2.3.1. Approximate input-state linearizing control This approach is an extension of the method presented in (Kanter et al., 2002); it involves requesting higher-order linear responses for the state variables (Panjapornpon et al., 2006). The following constrained optimization problem is solved at each time instant:  n

(τi D + 1)Pi xi − xssi 2 (7) min wi u τi Pi i=1

subject to: u(l) = 0,

l≥1

where xssi is the desired steady-state value of state, xi , corresponding to a given output set-point, ysp ; P1 ≥ R1 , . . ., Pn ≥ Rn ; R1 , . . ., Rn are relative orders of the state variables, x1 , . . ., xn , respectively; τ 1 , . . ., τ n are the time constants of the state variables in closed loop; w1 , ..., wn are adjustable, positive, scalar weights, whose values are set according to the relative importance of the state variables; as wi increases, the importance of the state variable xi increases. Using the model of (1), the time derivatives of the states in the performance index of (7) are again calculated analytically, which requires taking partial derivatives with respect to x and performing matrix operations. Note that the time derivatives of a state variable xi that are of order higher than Ri depend on the time derivative(s) of u. In the controller design, these time derivatives of u are set to zero, in principle, leading to a feedforward/static state feedback. 2.3.2. Input–output linearization with stability constraint This approach involves tracking linear output trajectories subject to a hard stability constraint at each time instant (Panjapornpon & Soroush, 2006): m 

(εi D + 1)ri yi − yspi 2 min (8) u εri i i=1

subject to: β2 V¨ + 2βV˙ + V ≤ 0

1574

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

Fig. 1. Flow diagram of the major components of the model-based controller design software.

where V is a Lyapunov function given by: V = (x − xss )T P(x − xss ) ε1 , . . ., εm are the time constants of the process outputs, y1 , . . ., ym , in closed loop, P is a positive-definite matrix, and the tuning parameter, β, is chosen such that β > max(ε1 , . . ., εm ). 3. Software environment DGMBC 1.0 has three principal components: MATHEMATICA routines, a simulation algorithm, and a user interface. A flow diagram of the software, showing the interactions among these components, is presented in Fig. 1. DGMBC 1.0 receives the process model and parameter values through the Front-End window, and then sends a set of commands to the MathLink program. The latter invokes the MATHEMATICA kernel to execute the desired control packages. The kernel handles the controller design task based upon the Process Information and the control method chosen by the user. The designed controller equations and user-chosen, tuning parameter values are used to carry out numerical simulations. The performance of the closed-loop control system is presented in the form of tables and plots. Each feature is described next. 3.1. MATHEMATICA packages DGMBC 1.0 has two types of MATHEMATICA packages, process analysis and controller design packages, to perform the controller design. The process analysis package named analyse.m consists of two main functions, StateEquilibrium and analyse. The StateEquilibrium function calculates feasible steady-state values within the given range of the state variables by solving the set of algebraic equations: 0 = f (x, dn , u),

ysp = h(x)

where dn is the nominal value of the measured disturbances, using symbolic inversion or Newton’s method with multiple initial guesses, as described in detail in Appendix A. The analyse function calculates the zero dynamics and examines the local asymptotic stability of the process and the process zero dynamics at the desired steady-state. The zero dynamics equations are obtained symbolically by setting the controlled outputs y1 , . . ., ym to their desired constant set-points and all time derivatives of the outputs to zero. They characterize the internal dynamics of the process when the controlled outputs are maintained at given constant set-points: dx dt ysp

= =

h(x)

0

=

h i (x)

0

=

h i (x, u(0) , u(1) , . . .)

f (x, dn , u) (9) i = 1, . . . , m, = 1, . . . , ri − 1 i = 1, . . . , m,

= ri , . . .

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

1575

For processes with a nonsingular characteristic matrix, the zero dynamics are of order (n − (r1 + · · · + rm )), the dimension of the space in which the system of (1) evolves when all controlled outputs are maintained at their constant set-points. The algorithm that MATHEMATICA uses to calculate the zero dynamics symbolically is described in Appendix A. The relative orders of the controlled outputs and the state variables are calculated using an algorithm based on the definition of the relative orders. For example, in the case of processes with no measured disturbances, to calculate the relative order of y1 , first ∂[dy1 /dt]/∂u is calculated symbolically:   ∂ dy1 ∂ ∂h1 (x) = f (x, u) ∂u dt ∂u ∂x If the right-hand side of the preceding equation is nonzero (if y˙ 1 depends on u), then r1 = 1. Otherwise, ∂[d2 y1 /dt2 ]/∂u is calculated symbolically: 

  ∂ d 2 y1 ∂ ∂ ∂h1 (x) = f (x, u) f (x, u) ∂u dt 2 ∂u ∂x ∂x If the right-hand side of the preceding equation is nonzero, then r1 = 2. Otherwise, this procedure is continued until the relative order is calculated. The controller design packages terminate the calculation of the relative orders when a relative order less than or equal the number of state variables cannot be found. The packages calculate the characteristic matrix symbolically after the relative orders have been calculated. The analyse function also evaluates the eigenvalues of the Jacobian of the zero dynamics and the Jacobian of the process at the given steady-state values. With the use of the MATHEMATICA kernel, the Jacobians of the process and the zero dynamics are calculated by symbolically differentiating: (a) the process state equations with respect to the process state variables, and (b) the zero dynamics equations with respect to the zero dynamics state variables, respectively. DGMBC 1.0 contains a collection of controller design packages that are based on the differential-geometric control methods presented in (Daoutidis et al., 1990; Kanter et al., 2002; Kravaris & Soroush, 1990; Panjapornpon & Soroush, 2006; Panjapornpon et al., 2006). Each of the controller design packages derives a controller consisting of a state feedback, a dynamic system to add integral action, and a state observer. It also generates the controller code in three programming languages: C, FORTRAN and MATLAB. To calculate of the state feedback, the time derivatives of the controlled outputs (state variables) up to the corresponding relative orders are calculated symbolically using the process model. These derivatives are substituted for in the relevant Eq. (2), (4), (7), or (8), and the equation is solved for u either symbolically or numerically. This solution involves symbolically or numerically solving a set of algebraic equations or numerically solving an optimization problem, as in (Panjapornpon & Soroush, 2006; Panjapornpon et al., 2006). The controller design packages such as Minimum.m, MinimumDm.m, and StableNonMinimum.m calculate first the relative orders of the outputs and then the state feedback analytically. The controller design packages such as UnstableNonminimum.m, and MIMOLyapunov.m formulate an optimization problem and its constraints, as described in (Panjapornpon & Soroush, 2006; Panjapornpon et al., 2006). The state feedbacks are then obtained by solving the corresponding optimization problem, (7) or (8), using the fmincon or fminunc function of MATLAB. The relevant optimization function is also integrated into the simulation algorithm described in the next section. 3.2. Simulation algorithm While MATHEMATICA has powerful symbolic manipulation capabilities, its numerical algorithms are less powerful. When a simulation algorithm retains huge data arrays and has many complicated decision loops, the computational speed of MATHEMATICA declines considerably. Because of this, the simulation algorithms of DGMBC were created as Dynamic Library Link (DLL) functions, creating a flexible environment for a stand-alone application. The DLL files developed by the Component Object Model (COM) Builder of MATLAB contains a plot function, an ODE solver function, ode45, and two optimization functions, fminunc and fmincon. The controller equations in the MATLAB format, with given process parameters, are integrated by the simulation functions, OdeAna, OdeFminunc and OdeFmincon, that were developed specifically for the control methods in (Daoutidis et al., 1990; Kanter et al., 2002; Kravaris & Soroush, 1990; Panjapornpon & Soroush, 2006; Panjapornpon et al., 2006) under a perfect-model assumption. While DGMBC 1.0 is for unconstrained processes, the next version of the software will be applicable to constrained processes and will allow for multiple set-point changes, step changes in parameter values, addition of measurement noise, and instantaneous perturbations in the state values. 3.3. User interface The Front-End window interacts with the user in five steps: Process Information, Process Analyses, Controller Equations, Closedloop Simulation, and Plot Selection. The user enters specific information in each step before proceeding to the next step. A flow diagram linking the steps and the tasks that DGMBC 1.0 performs in each step is shown in Fig. 2. In addition, the Front-End window

1576

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

Fig. 2. Flow diagram of the tasks in the model-based controller design software (numbers correspond to the steps in the user interface).

has the Input Command step that shows all of the commands sent to the MATHEMATICA kernel. The user enters the process model equations and identifies the controlled outputs, manipulated inputs, process parameters, and available measurements in the Process Information step shown in Fig. 3. The process model is entered as a set of ordinary differential equations in the MATHEMATICA format. In the Process Analyses step (Fig. 4a), the software calculates feasible steady-state values within the given variable ranges. The user then selects a desired steady-state pair from those calculated by the software, as illustrated in Fig. 4b. Upon clicking on the Calculate button, DGMBC 1.0 calculates the zero dynamics of the process, the eigenvalues of the Jacobian of the process evaluated at the desired steady-state, and the eigenvalues of the Jacobian of the process zero dynamics evaluated at the desired steady-state. Depending on the analysis results (stability and minimum-phase-ness), DGMBC 1.0 allows the user to select an appropriate controller method. According to the logic shown in Fig. 5, the software directs the user to suitable control methods. In the case that the user selects: (a) approximate input–output linearization, the order of the response of each controlled output should be chosen; and (b) approximate input-state linearization, the order of the response of each state variable should be chosen. The CodeViewer window in the Controller Equations step (Fig. 6) presents the model-based controller equations in the user’s desired formats, C, FORTRAN and MATLAB (Fig. 7). Subsequently, a file containing the code for the control (closed-loop) system is loaded in the Closed-Loop Simulation step (Fig. 8). The user enters initial values for the state variables, controller-tuning parameter values, and the simulation time. The simulation method in the dynamic library link (DLL) is selected according to the type of the control method. The DDL is a library of executable functions that can be used by a Windows application. The creation of the simulation functions as DLLs allows one to install the new software package on any Windows-based computer without requiring the installation of MATLAB. The closed-loop simulation results are presented in the form of tables and plots. The user selects the plot variables in the Plot Selection step (Fig. 9). The user can save the entered information in the Process Information, Process Analyses and Closed-loop Simulation steps by clicking on Save or Save As (in the File menu), which saves the information in each of the steps in a file in “Model-Based Input” format (*.mi). He/she can also load the saved information back to DGMBC 1.0 by clicking on Open in the File menu. 4. Examples This section presents the application of the software to three process examples.

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

1577

Fig. 3. The Process Information step with the bioreactor example information.

4.1. Continuous bioreactor Consider a constant-volume, continuous bioreactor described by:   x2 x˙ 1 = −x1 u + x1 (1 − x2 ) exp γ  1 1 + α1 x2 x˙ 2 = −x2 u + x1 (1 − x2 ) exp γ1 1 + α 1 − x 2 y = x1 where x1 is the dimensionless outlet cell mass concentration, x2 is the dimensionless outlet substrate concentration, and u is the dimensionless substrate feed rate. The parameter values are as follows: α1 = 0.02 and γ 1 = 0.48. It is desired to maintain the cell mass concentration, x1 , at 0.1448 by manipulating u. It is assumed that only the cell mass concentration is measured. With x = [x1 x2 ]T , u = [u], and y = [x1 ], as shown in Fig. 3, the following process information is entered in the Process Information step: 4.1.1. Process information Process model x1 [t ] = −x1[t] ∗ u[t] + x1[t] ∗ (1 − x2[t]) ∗ Exp[x2[t]/gamma1] x2 [t ] = −x2[t] ∗ u[t] + x1[t] ∗ (1 − x2[t]) ∗ Exp[x2[t]/gamma1] ∗ (1 + alpha1)/(1 + alpha1 − x2[t]) y1 = x1[t] Process parameter value(s): alpha1 = 0.02, gamma1 = 0.48 State variable(s): x1[t], x2[t] Controlled output(s): y1 State and manipulated input ranges: {x1[t],0,1},{x2[t],0,1},{u[t], 0.2, 1.5} Interval of search ranges: 5 Option: Observer design Option: Continuous process

Manipulated input(s): u[t] Output setpoint value(s): 0.1448

1578

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

Fig. 4. (a) The Process Analyses step with the bioreactor example information. (b) The calculated feasible steady-state for the bioreactor.

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

1579

Fig. 5. Four classes (quadrants) of processes and suitable controller methods for each quadrant. IOL = Input-output linearization (Subsection 2.1); AIOL = approximate input–output linearization (Subsection 2.2); AISL = approximate input-state linearization (Subsection 2.3.1); IOLSC = input-output linearization with a stability constraint (Subsection 2.3.2).

Upon clicking on the right arrow button in the bottom left of the Process Information window, DGMBC 1.0 calculates feasible steady-state pairs within the given operating ranges and presents the results in the Process Analyses step, shown in Fig. 4a. The user then chooses and enters the pair (x1ss = 0.1448, x2ss = 0.8455, uss = 0.9) as the desired steady-state pair. Upon clicking on the Calculate button (Fig. 4a), the software calculates the zero dynamics of the process, the eigenvalues of the Jacobian of the process evaluated at the desired steady-state, and the eigenvalues of the Jacobian of the process zero dynamics evaluated at the desired steady-state: λ1 (Jol ) = 0.061 + 1.731i, λ2 (Jol ) = 0.061 − 1.731i and λ1 (Jzd ) =3.449 (Fig. 4b). Because the desired steady-state pair is unstable and non-minimum-phase, DGMBC 1.0 offers two choices for the controller design method, approximate input–output linearization (Panjapornpon et al., 2006) and input–output linearization with stability constraints (Panjapornpon & Soroush, 2006).

Fig. 6. The generated controller equations for the bioreactor in Controller Equations step.

1580

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

Fig. 7. The CodeViewer window showing controller equations in the C language format.

Fig. 8. Closed-Loop Simulation step for the bioreactor.

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

1581

Fig. 9. The Plot Selection step for the bioreactor.

The approximate input-state linearizing control method (Panjapornpon et al., 2006) is selected as the desired controller design method, and the orders of the requested state responses are set: p1 = 1 and p2 = 1 (at the bottom of the Process Analysis window). The information entered in the Process Analyses step is given below and in Fig. 4b. 4.1.2. Process analyses Desired steady-state pair: 0.1448, 0.8455, 0.9 Controller methods: Approximate input–output linearization

Order(s) of output response(s): 1, 1

Upon clicking on the right arrow at the bottom of the window (Fig. 4b), the software generates the controller equations shown in Fig. 6. The Controller Equations step allows one to obtain the controller equations in the various formats. After choosing (clicking on) the desired code format at the bottom of the window and then clicking on the View button, DGMBC 1.0 shows the controller code in the CodeViewer window (Fig. 7). Clicking on the right arrow at the bottom of the window of Controller Equations step (Fig. 6) takes the user to the Closed-loop Simulation step (Fig. 8). In the Closed-Loop Simulation step, the initial conditions and the controller-tunable-parameter values are entered: [x1 (0), x2 (0), ζ1 (0), w1 (0), w2 (0)] = [0.1, 0.75, 0.05, 0.1, 0.75], ε1 = 0.1, ε2 = 0.35, and L2 = 0.7. A simulation time of 4 h and an initial guess of 0.45 for the manipulated variable are chosen. 4.1.3. Closed-loop simulation Simulation time: [0, 4] Tuning parameter(s): {Epsilon1, Epsilon2, L11} State variable(s) of the closed-loop system: x1, x2, Zeta1, w1, w2 Initial condition(s): 0.1, 0.75, 0.05, 0.1, 0.75 Manipulated variable(s): u

Tuning value(s): 0.1, 0.35, 0.7

Initial guess for the manipulated variable(s): 0.45

Upon clicking on the Simulation button located in the top right corner of the window, DGMBC simulates the closed-loop system. Clicking on the right arrow button at the bottom of the Closed-loop Simulation window takes the user to the Plot Selection step (window). The axis choices shown in the Plot Selection step of Fig. 9 are made using the dropdown menus. After clicking the Plot button located at the bottom right corner of the window, the Fig. 1 window shown in Fig. 10 opens. The figure shows the closed-loop responses of the variables selected in the Plot Selection step.

1582

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

Fig. 10. Simulated closed-loop responses of the controlled output and manipulated input of the bioreactor.

4.2. Multi-input, multi-output jacketed chemical reactor Consider a non-isothermal, continuous chemical reactor in which the reactions A → B → C take place in the liquid phase. The process dynamics are represented by the following model: dCA dt dCB dt dT dt dTj dt

  E1 F 2 = (CA0 − CA ) − k1 exp − CA V RT     F E1 E2 2 = − CB + k1 exp − − k2 exp − CA CB V RT RT     (−H1 ) E1 (−H2 ) E2 US F 2 k1 exp − + k2 exp − = (Ti − T ) + CA CB + (Tj − T ) V ρcp RT ρcp RT ρcp V Fj US = (Tji − Tj ) − (Tj − T ) V ρj cpj Vj

y1 = CB y2 = T It is desired to maintain CB at 5.233 mol/l and T at 443.92 K by manipulating F and Fj . Only measurements of CB and T are available. Let x = [CA CB T Tj ]T , u = [F Fj ]T and y = [CB T]T . The process dynamics are described by:   u1 E1 x12 (CA0 − x1 ) − k1 exp − V Rx3     u1 E1 E2 x12 − k2 exp − x2 x˙ 2 = − x2 + k1 exp − V Rx3 Rx3     u1 (−H1 ) E1 (−H2 ) E2 US 2 x1 + x2 + x˙ 3 = k1 exp − k2 exp − (Ti − x3 ) + (x4 − x3 ) V ρcp Rx3 ρcp Rx3 ρcp V u2 US x˙ 4 = (Tji − x4 ) − (x4 − x3 ) V ρj cpj Vj x˙ 1 =

y1 = x2 y2 = x3 The information entered in the Process Information step of DGMBC 1.0 is given below.

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

1583

4.2.1. Process information Process model x1 [t ] = −((k1 ∗ x1[t]ˆ2)/Eˆ(E1R/x3[t])) + ((Ca0 − x1[t]) ∗ u1[t])/V x2 [t ] = (k1 ∗ x1[t]ˆ2)/Eˆ(E1R/x3[t]) − (k2 ∗ x2[t])/Eˆ(E2R/x3[t]) − (x2[t] ∗ u1[t])/V x3 [t ] = −((k1 ∗ H1 ∗ x1[t]ˆ2)/(Cp ∗ Eˆ(E1R/x3[t]) ∗ Rho)) − (k2 ∗ H2 ∗ x2[t])/(Cp ∗ Eˆ(E2R/x3[t]) ∗ Rho) + (S ∗ U ∗ (−x3[t] + x4[t]))/(Cp ∗ Rho ∗ V) + ((Ti − x3[t]) ∗ u1[t])/V  x4 [t ] = −((S ∗ U ∗ (−x3[t] + x4[t]))/(Cpj ∗ Rhoj ∗ Vj)) + ((TJi − x4[t]) ∗ u2[t])/Vj y1 = x2[t] y2 = x3[t] Process parameter value(s): k1 = 2.5*10 ˆ 10, k2 = 1.5*10 ˆ 10, H1 = −20, H2 = −80, Ca0 = 12, V = 5, Vj = 5, E1R = 8000, E2R = 9100, Rho = 1, Rhoj = 1.1, Cp = 2.25, Cpj = 3, Ti = 320, TJi = 298.15, U = 3825, S = 0.225 State variable(s): x1[t], x2[t], x3[t], x4[t] Manipulated input(s): u1[t], u2[t] Controlled output(s): y1, y2 Output setpoint value(s): 5.233, 443.92 State and manipulated input range: {x1[t],0,1.5}, {x2[t],0,10}, {x3[t],400,450}, {x4[t],400,440}, {u1[t],50,120}, {u2[t],50,120} Interval of search ranges: 5 Option: Observer design Option: Continuous process

The steady-state pair corresponding to the desired set-point, y1sp = 5.233, and y2sp = 443.92 (x1ss = 0.7, x2ss = 5.233, x3ss = 443.92, x4ss =403.24, u1ss = 80.95, and u2ss =100.94), is selected. The process is unstable [λ1 (Jol ) = −491.38, λ2 (Jol ) = −90.71, λ3 (Jol ) =86.28, λ4 (Jol ) =−15.36] and non-minimum-phase [λ1 (Jzd ) = 589.59]. The approximate input-state linearizing control method (Panjapornpon et al., 2006) is selected as the controller design method. The orders of the linear responses of state variables, x1 , x2 , x3 , and x4 , are requested to be p1 = 2, p2 = 2, p3 = 2, and p4 = 1, respectively. The inputs in the Process Analyses step are given below. 4.2.2. Process analyses Desired steady-state pair: 0.700, 5.233, 443.92, 403.24, 80.95, 100.94 Order(s) of output response(s): 2, 2, 2, 1 Controller methods: Approximate input-state linearization

The process is initially at [x1 (0), x2 (0), x3 (0), x4 (0)] = [0.7, 5.233, 443.92, 403.24], y1sp = 5.233, and y2sp =443.92. The set-point is in the non-minimum-phase region. The following values are selected: ε1 = 0.09, ε2 = 0.09, ε3 = 0.08, ε4 = 0.08, L1 = 0.5, L2 = 0.5, with a simulation time of 3 h. The inputs in the Closed-Loop Simulation step are given below. The simulation results showing the controller performance are shown in Fig. 11. 4.2.3. Closed-loop simulation Simulation time: [0,3] Tuning parameter(s): Epsilon1, Epsilon2, Epsilon3, Epsilon4, L11, L22 Tuning value(s): 0.09, 0.09, 0.08, 0.08, 0.5,0.5 State variable(s) of the closed-loop system: x1, x2, x3, x4, Zeta1, Zeta2, w1, w2, w3, w4 Initial condition(s): 1.8138, 9.4456, 387.69, 380.67, −2.909, 186.83, 1.8138, −5.8085, 9.4456, 2.8518, 387.69, 25.457, 380.67 Manipulated variable(s): u1, u2

Initial guess for the manipulated variable(s): 40, 10

4.3. Continuous chemical reactors in series Consider two isothermal, continuous stirred tank reactors in series, in which the reactions A → B → C take place in the liquid phase. The process model is: dCA1 dt dCB1 dt dCA2 dt dCB2 dt y1

F1 2 (CA0 − CA1 ) − kA CA1 V1 F1 2 − kB CB1 = − CB1 + kA CA1 V1 F2 F1 2 (CA1 − CA2 ) + (CA0 − CA2 ) − kA CA2 = V2 V2 F2 F1 2 (CB1 − CB2 ) − CB2 + kA CA2 − kB CB2 = V2 V2 = CB1

y2

= CB2

=

1584

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

Fig. 11. Simulated closed-loop responses of the controlled outputs and manipulated inputs of the MIMO chemical reactor.

where F1 and F2 are the inlet volumetric flow rates of pure A, V1 = 0.01, V2 = 0.01, CA0 = 7, kA = 6, and kB = 1. The concentrations, CB1 and CB2 , are measured. It is desired to maintain these two concentrations at y1sp = 2, and y2sp = 3 by manipulating the feed rates, F1 and F2 . The reactor parameter values are given in Table 1. With x = [CA1 CB1 CA2 CB2 ]T , u = [F1 F2 ]T , and y = [CB1 CB2 ]T , the inputs in the Process Information step are given below.

4.3.1. Process information Process model x1 [t ] = −(ka ∗ x1[t]ˆ2) + ((Ca0 − x1[t]) ∗ u1[t])/V1 x2 [t ] = ka ∗ x1[t]ˆ2 − kb ∗ x2[t] − (x2[t] ∗ u1[t])/V1 x3 [t ] = −(ka ∗ x3[t]ˆ2) + ((x1[t] − x3[t]) ∗ u1[t])/V2 + ((Ca0 − x3[t]) ∗ u2[t])/V2 x4 [t ] = ka ∗ x3[t]ˆ2 − kb ∗ x4[t] + ((x2[t] − x4[t]) ∗ u1[t])/V2 − (x4[t] ∗ u2[t])/V2 y1 = x2[t] y2 = x4[t]

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

1585

Table 1 Process parameters of the non-isothermal jacketed reactor Parameter

Value

Unit

CA0 k1 k2 E1 /R E2 /R Ti Tji H1 H2 ρ ρj cp cpj U S V Vj

12 2.5 × 1010 1.5 × 1010 8,000 9,100 320 298.15 −20 −80 1 1.1 2.25 3 3,825 0.225 5 5

mol l−1 l mol−1 h−1 h−1 K K K K kJ mol−1 kJ mol−1 kg l−1 kg l−1 kJ kg−1 K−1 kJ kg−1 K−1 kJ m−2 K−1 s−1 m2 l l

Fig. 12. Simulated closed-loop responses of the controlled outputs and manipulated inputs of the isothermal chemical reactors in series.

1586

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

Process parameter value(s): ka = 6, kb = 1, Ca0 = 7, V1 = 0.01, V2 = 0.01 State variable(s): x1[t], x2[t], x3[t], x4[t] Controlled output(s): y1, y2 State and manipulated input range: {x1[t],0,7},{x2[t],0,7},{x3[t],0,7},{x4[t],0,7}, {u1[t],0,0.2}, {u2[t],0,0.2} Interval of search ranges: 5 Option: Observer design Option: Continuous process

Manipulated input(s): u1[t],u2[t] Output setpoint value(s): 2, 3

The steady-state pair (x1ss = 0.699, x2ss = 2, x3ss = 1.101, x4ss = 3, u1ss = 0.005, u2ss = 0.0013) that corresponds to y1sp = 2 and y2sp = 3 is selected. The pair is stable [λ1 (Jol ) = −14.90, λ2 (Jol ) = −8.85, λ3 (Jol ) = −2.72, λ4 (Jol ) = −1.46] and non-minimum-phase [λ1 (Jzd ) = 17.57, λ2 (Jzd ) = 11.03]. The approximate input–output linearizing control method (Kanter et al., 2002), with p1 = 2 and p2 = 2, is used. The inputs in the Process Analyses step are detailed below. 4.3.2. Process analyses Desired steady-state pair: 0.699, 2.0, 1.101, 3.0, 0.005, 0.013 Order(s) of output response(s): 2, 2 Controller methods: Approximate input–output linearization

With [x1 (0), x2 (0), x3 (0), x4 (0)] = [0.75, 2.192, 0.75, 2.192], ε1 = 0.4, ε2 = 0.3, and 6 h of simulation time, the inputs in the Closed-loop Simulation step are given below. The closed-loop responses are shown in Fig. 12. 4.3.3. Closed-loop simulation Simulation time: [0,6] Tuning parameter(s): {Epsilon1, Epsilon2} State variable(s) of the closed-loop system: {x1, x2, x3, x4, w1,w2,w3,w4} Initial condition(s): 0.75,2.192,0.75,2.192,0.75,2.192,0.75,2.192 Manipulated variable(s): {u1, u2}

Tuning value(s): 0.4, 0.3

Initial guess for the manipulated variable(s): 0.1, 0.1

5. Conclusions DGMBC 1.0, a user-friendly, integrated, software package for the design of differential-geometric model-based controllers, was presented. Given a lumped-parameter process model, DGMBC derives an analytical model-based controller, generates the controller equations in C, FORTRAN and MATLAB format, and carries out closed-loop simulations. DGMBC allows control engineers to design differential-geometric, model-based controllers with ease. With the software, differential-geometric, model-based controllers can be designed for processes involving complex models, as implemented in industry. A MathLink program provides an interface for communication with the MATHEMATICA kernel. DGMBC 1.0 was applied to design model-based controllers for chemical and biochemical reactors. A soft copy of DGMBC 1.0 can be obtained by contacting the corresponding author of this paper. The next version of the software will be applicable to constrained processes and will allow for multiple set-point changes, steps in parameter values, and instantaneous perturbations in the state values. Acknowledgements This study was supported by the National Science Foundation Grants CTS-0101133 and CTS-0101237. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. Appendix A A.1. StateEquilibrium function Under the following assumptions made in Section 2: (a) The Jacobian of the process, ∂f(x,d,u)/∂x, and [∂h(x)/∂x][∂f(x,d,u)/∂x]−1 [∂f(x,d,u)/∂u] evaluated at the desired steady-state are nonsingular. (b) The numbers of inputs and outputs are equal. 1. Set the process parameters to their values in the equations. 2. Solve symbolically the set of equations f (x, dn , u) = 0,

ysp = h(x)

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

1587

for x and u by using the Solve function of MATHEMATICA. If the equations can be solved by the Solve function (no error message from MATHEMATICA), then go to step 4. If not, the StateEquilibrium function then switches to the FindRoot function of MATHEMATICA (step 3), which uses Newton’s method. 3. Newton’s Method (a) Discretize the range of each variable into eight equal intervals (eight is the default value, but the user can change this number): Given xi ∈ [xil , xih ] and uj ∈ [ujl , ujh ], divide the xi range into the 8 equal intervals [xil , xil + Δxi ], . . . , [xih − Δxi , xih ] and the uj range into the 8 equal intervals [ujl , ujl + Δuj ], . . . , [ujh − Δuj , ujh ], where Δxi =

xih − xil uj − ujl and Δuj = h 8 8

(b) In each of the 8nm hyper-cubes formed by the intervals, select an initial guess (xguess , uguess ) for Newton’s method. Run the FindRoot function with each of the initial guesses, totally 8 nm times. (c) Screen all solutions from the 8 nm times in step (b) to obtain the feasible equilibrium pairs (those that satisfy xi ∈ [xil , xih ] and uj ∈ [ujl , ujh ], ∀i, j) and then go to step 4. 4. Show the solution in the user interface. A.2. Analyse function With the same assumptions made for the StateEquilibrium function in the preceding section, given the desired steady-state (xss , dss , uss ), set the process parameters to their values in the equations. 1. Check the asymptotic stability of the process: (a) Symbolically calculate the Jacobian of the process using MATHEMATICA (b) Calculate the eigenvalues of the Jacobian evaluated at the desired steady-state using the Eigenvalues function of MATHEMATICA 2. Calculate the zero dynamics: (a) Calculate the relative order of each controlled variable and the characteristic matrix using MATHEMATICA (i) Calculate the ri th derivative of yi (i = 1, . . ., m) symbolically using MATHEMATICA (ii) Calculate the characteristic matrix of the process symbolically using MATHEMATICA (b) Calculate the unforced zero dynamics. The zero dynamics describe the dynamics of the process when all controlled outputs are maintained at constant set-points. In the case of processes without measured disturbances (d = 0), the following set of algebraic equations hold (Isidori, 1989): yspi = hi (x), i = 1, . . . , m (j)

0 = yi :=hj (x), 0 = yi(ri ) = hri (x, u),

i = 1, . . . , m, j = 1, . . . , (ri − 1) i = 1, . . . , m,

( ) yi

0= = h (x, u, u(1) , ...), i = 1, . . . , m, = ri + 1, . . . When the characteristic matrix is nonsingular, the preceding set of algebraic equations is adequate for calculating the zero dynamics (Isidori, 1989): yspi = hi (x), i = 1, . . . , m, (j)

0 = yi :=hj (x), i = 1, . . . , m, j = 1, . . . , (ri − 1) (r ) r 0 = yi i = h i (x, u), i = 1, . . . , m For processes with a nonsingular characteristic matrix, MATHEMATICA solves the preceding system of algebraic equations symbolically and obtains the manipulated inputs and (r1 + · · · + rm ) of the state variables in terms of the remaining (n − (r1 + · · · + rm )) state variables. Then, it substitutes for the (r1 + · · · + rm ) state variables and the manipulated inputs in the state equations of the model of (1) to obtain the dynamic system describing the (n − (r1 + · · · + rm )) state variables. An algorithm for calculating the zero dynamics of processes with a singular characteristic matrix is presented in (Isidori, 1989). References Barker, H. A., & Zhuang, M. (1997). Control system analysis using Mathematica and a graphical user interface. Computing & Control Engineering Journal, 8, 64–69. Barker, H. A., Grant, P. W., & Zhuang, M. (1996). Control system analysis with Mathematica. IEE Colloquium on Symbolic Computation for Control, 2, 1–5. Bemporad, A., & Morari, M. (2001). Optimization-based hybrid control tools. In Proceedings of American Control Conference, vol. 2 (pp. 1689–1703). Blankenship, G. L., Kwatny, H. G., LaVigna, C., & Polyakov, V. (1997). Integrated modeling and design of nonlinear control systems. In Proceedings of American Control Conference, vol. 3 (pp. 1395–1399).

1588

C. Panjapornpon et al. / Computers and Chemical Engineering 32 (2008) 1569–1588

Campbell, S. L., Delebecque, F., & von Wissel, D. (1996). A mixed symbolic-numeric software environment. In Proceedings of IEEE International Symposium on Computer-Aided Control System Design (pp. 436–441). Daoutidis, P., Soroush, M., & Kravaris, C. (1990). Feedforward/feedback control of multivariable nonlinear processes. American Institute of Chemical Engineers Journal, 36, 1471–1484. Datta, B. N., & Sarkissian, D. (1999). Numerical linear control library; a Mathematica-based integrated environment with the modern control algorithms. In Proceedings of IEEE International Symposium on Computer-Aided Control System Design (pp. 91–96). de Jager, B. (1995). The use of symbolic computation in nonlinear control: Is it viable? IEEE Transactions on Automatic Control, 40, 84–89. Isidori, A. (1989). Nonlinear control systems: An introduction (2nd ed.). Berlin: Springer-Verlag. Kanter, J. M., Soroush, M., & Seider, W. D. (2002). Nonlinear feedback control of multivariable non-minimum-phase processes. Journal of Process Control, 12, 667–686. Kitamoto, T., Saeki, M., & Ando, K. (1992). CAD package for control system on Mathematica. In IEEE International Conference on Systems Engineering (pp. 448–451). Kravaris, C., & Soroush, M. (1990). Synthesis of multivariable nonlinear controllers by input/output linearization. American Institute of Chemical Engineers Journal, 36, 249–264. Maple6. (2000). Waterloo. Canada: Waterloo Maple Inc. Mathematica. (1999). MATHEMATICA (4th ed.). Champaign, IL: Wolfram Research Inc. Matlab. (2002). MATLAB., 7. 0 Ed. Natick, MA: The Mathworks Inc. Meel, A., Seider, W. D., & Soroush, M. (2006). A Game theoretic approach to multi-objective process designs: Focus on inherent safety. American Institute of Chemical Engineers Journal, 52(1), 228–246. Ohtani, T., Fukuzawa, M., & Masubuchi, M. (1994). A CAD system for nonlinear control system design using Mathematica. In Proceedings of IEEE/IFAC Joint Symposium on Computer-Aided Control System Design (pp. 197–204). Panjapornpon, C., & Soroush, M. (2006). Control of non-minimum-phase nonlinear systems through constrained input-output linearization. In Proceedings of American Control Conference (pp. 4522–4527). Panjapornpon, C., Soroush, M., & Seider, W. D. (2006). Model-based controller design for unstable, non-minimum-phase, nonlinear processes. Industrial and Engineering Chemistry Research, 45, 2758–2768. Sack, J., & Singh, T. (2000). Automated design of model predictive controllers. In Proceedings of American Control Conference (pp. 3758–3762). Soroush, M. (1998). State and parameter estimations and their applications in process control. Computer and Chemical Engineering, 23, 229–245.