Optimal control of hybrid dynamical systems: application in process engineering

Control Engineering Practice 10 (2002) 133–149

Optimal control of hybrid dynamical systems: application in process engineering Philippe Manon, Claire Valentin-Roubinet*, Ge! rard Gilles # 308G, CPE, 43 Bd du 11 Novembre 1918, 69622 Villeurbanne, France LAGEP, Universite! Claude Bernard Lyon1, Bat. Received 27 July 2000; accepted 8 June 2001

Abstract This paper deals with a class of hybrid dynamical systems which may operate in multiple phases executed in parallel or in sequence. Within each phase the process dynamic is continuous. Two different methods of designing a control system which minimizes the overall operating time used in driving the process from a default mode to one of the nominal modes, respecting the constraints on the continuous variables, are presented. The first is based on minimization of a criterion under constraints. The other is based on Pontryagin’s maximum principle extended to discrete event controlled hybrid dynamical systems. Both are illustrated on a product supply process. r 2002 Published by Elsevier Science Ltd. Keywords: Hybrid dynamical systems; Hybrid automaton; Non-linear optimal control; Pontryagin’s maximum principle; Optimization; Chemical process

1. Introduction Chemical processes include manufacturing phases such as, filling or emptying a reactor, heating or mixing a product. They may operate in multiple distinct phases (or modes) which can be executed in parallel or in sequence. The different phases need several resources such as valves, pipes, reactors, heaters, which have characteristics of different natures: actuators admit a logic control (to open, to close, to heat, to cool, etc.) and sensors give a continuous information (flow rate, heating power, temperature, etc.). The model for such processes is then hybrid in the sense that it includes both discrete and continuous variables. The sequences of phases are represented by discrete or Boolean variables, and the dynamic evolution of levels, concentrations or temperatures in the equipment is represented by differential and/or algebraic equations. The structure of the controlled hybrid system is shown in Fig. 1. This paper presents two different methods of designing a control system which minimizes the overall time

*Corresponding author. Fax: +33-4-7243-1699. E-mail address: [email protected], [email protected] (C. Valentin-Roubinet). 0967-0661/02/$ - see front matter r 2002 Published by Elsevier Science Ltd. PII: S 0 9 6 7 - 0 6 6 1 ( 0 1 ) 0 0 1 2 3 - X

duration necessary for the system to reach a final hybrid region (a hybrid region is made of a region in the continuous state space and a discrete phase) from an initial hybrid region, respecting the constraints on the continuous variables and the allowed discrete event control sequences. If the controllable sequence to be optimized is not identified, it can be systematically determined by a controllability analysis method (Manon & Valentin-Roubinet, 1999a, b). In this paper, only the sequence which has a minimum number of phases is optimized. The process is modelled by a hybrid automaton with a discrete event control. The control hybrid sequence is composed of a sequence of phases (or discrete states), each with its own time duration. The sequential controller (cf. Fig. 1) generates this hybrid sequence which consists of discrete controllable events impulsed at precisely precalculated time instants. The first method explains how an optimal control problem for hybrid dynamical systems (HDS) can be translated in terms of a functional and constraints set. The second method shows how Pontryagin’s maximum principle (PMP) can be extended to HDS with a discrete event control. Section 2 presents the hybrid automaton model and is applied to a product supply process.

134

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

Fig. 1. Closed loop control structure of the hybrid dynamic system.

Section 3 presents the first optimization method. First the problem position is made clearer and the existence conditions of an optimal solution is investigated. Then, the translation in terms of functional and constraints are presented and this method is applied to the product supply process. Section 4 presents the second optimal control synthesis method. First, it is verified that the PMP can be applied to HDS with a discrete event control. Next, the necessary conditions for a control to be optimal are given. These conditions lead to a differential-algebraic system of equations to be solved. Then, the method is applied to the same process as presented in Section 3. Finally, the two methods are compared.

denoted by l1 ¼ fl1 ðlj Þ; jAf1; y; mgg; that is 8lj AQ; x’ ¼ l1 ðlj ÞðxÞ:The function l1 ðlj Þ is supposed to be locally Lipschitz with respect to x: * m 2 is the set of constraints for each phase: m2 ¼ fm2 ðlj Þ=jAf1; y; mgg: m2 ðlj Þ represents the constraints for the phase lj . It is described by cj linear inequalities: 8kAf1; y; cj g; CkT xpdk with Ck an nlength vector and dk a constant. It is not necessarily a state-space partition. The constraints m2 ðlj Þ cannot be violated, which means that phase lj is left as soon as an inequality of m2 ðlj Þ is going to be false. * S is the events set. S ¼ fs ; jAf1; y; mg; kA jk f1; y; mg=jakg: Each event sjk defines the discrete transition trjk AQ2 from location lj to location lk : A discrete transition trjk can be either autonomous or controlled. If it happens on constraints m2 ðlj Þ boundary, the event sjk is generated by the plant, sjk is uncontrolled and the transition is autonomous. Or, if it happens inside m2ðlj Þ region, sjk is generated by the control system, sjk is controlled and the transition is controlled. * m is a set of functions associated with the events s : jk 3 ( ) m3 ðlj ; sjk ; lk Þ= m3 ¼ jAf1; y; mg; kAf1; y; mg; jak; sjk is defined It defines real valued variable jumps when an event occurs: * if s jk exists, then the phase transition from lj to lk is defined and m3 ðlj ; sjk ; lk Þ : X-X=xðtþ Þ ¼ m3 ðlj ; sjk ; lk Þðxðt ÞÞ

2. Modelling 2.1. Hybrid modelling The hybrid process is modelled by a hybrid automaton (Tittus & Egardt, 1998; Asarin, Maler, & Pnueli, 1995). As the Branicky unified model (Branicky, Borkar, & Mitter, 1994), autonomous and controlled jumps (continuous variable discontinuity) and autonomous and controlled switching (vector field discontinuity) may be represented. Definition 1. A hybrid automaton HA is defined by an 8-tuple HA ¼ ðX; Q; m1 ; m2 ; S; m3 ; Q0; Qf Þ: *

*

*

X is a vector subspace XDRn which represents the state-space of the system. x is the continuous state vector of the system denoted by x ¼ ½x1 yxn T AX: It is supposed to be continuously observable. Q is the possible discrete phases set (also called modes, locations or discrete states) of the hybrid process: Q ¼ flj ; jAf1; :::; mgg: Then, the hybrid state of the system is given by the couple ðx; lÞAX  Q: l1 is the set of m vector fields associated with each phase lj : The dynamics of the state vector x are

if sjk does not exist, neither the phase transition from lj to lk ; nor m3 ðlj ; sjk ; lk Þ are defined. If m3 ðlj ; sjk ; lk Þðxðt ÞÞaxðt Þ; a state jump (or state discontinuity) occurs. Q0 is the set of initial locations and Qf the set of final locations. *

*

2.2. Modelling of a product supply process So as to illustrate the method, the product supply process sketched in Fig. 2 will be modelled. It is a continuously stirred tank reactor (i.e. the temperature and the concentration are homogeneous inside the reactor). The aim of the system is to offer the customer a product with the following given characteristics: * *

*

a given flow rate F; a concentration C remaining between minimal and maximal values, and a temperature T remaining between minimal and maximal values.

The incoming products are provided by two different suppliers by either the valve V1 or the valve V2 : Therefore, the incoming products have different con-

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

135

Fig. 3. Process model.

Fig. 2. Products supply process.

centrations (C1 or C2 ) and different temperatures (T1 or T2 ). When the reactor is filled via valve V1 ; a heating system at temperature Tj1 (defining the first dynamic phase) is coupled. When valve V2 is open, the heating system delivers a temperature Tj2 (defining the second dynamic phase). The heating system is realized by a jacket (temperature index j) which surrounds the reactor, a tank and a pump. The reactor characteristics are given below: * * *

3

maximum capacity: VM ¼ 10 m , thermal exchange area: A ¼ 28:2 m2, thermal exchange coefficient: U ¼ 72:09 J (s m2 K)1.

The liquid supplied at the output has the following properties: * * *

molecular weight: wm ¼ 0:01802 kg mol1, density: d ¼ 1000 kg m3, heat capacity: Cp ¼ 33:7 J. (K mol)1.

In the first dynamic phase, the incoming liquid has the following properties: * * * *

concentration: C1 ¼ 100 mol m3, flow rate: F1 ¼ 5  104 m3 s1, temperature: T1 ¼ 308 K, Temperature in the jacket: Tj1 ¼ 308 K.

The continuous state vector has 3 components: x1 is the liquid output concentration C; x2 is the liquid output temperature T and x3 is the liquid volume V inside the reactor. The differential non-linear equations in the dynamic phases, x’ ¼ l1 ðli ÞðxÞ where iAf1; 2; 3g (Corriou, 1996), are as follows: x’ 1 ¼

Fi ðCi  x1 Þ; x3

ð1Þ

x’ 2 ¼

Fi UAwm ðTi  x2 Þ þ ðTji  x2 Þ; x3 dCp x3

ð2Þ

x’ 3 ¼ Fi  F:

The system must be driven from an initial phase l0 ; where the continuous region is reduced to a single point x0 (m2 ðl0 Þ : x1 ¼ 2500 mol m3, x2 ¼ 303 K, x3 ¼ 10 m3) and the dynamics are null (i.e. x’ ¼ l1 ðl0 Þ ¼ 0), to the desired final phase (where dynamics are also null: x’ ¼ l1 ðl3 Þ ¼ 0) with the constraints m2 ðl3 Þ: 2200p x1 p2300; 311px2 p315; 0px3 p10: The global constraints GC are: 2000px1 p2500; 303px2 p323; 0px3 p10 and m2 ðl1 Þ ¼ m2 ðl2 Þ ¼ GC: The process model is given in Fig. 3. In order to simplify the hybrid automaton graphical representation, *

*

In the second dynamic phase, the incoming liquid has the following properties: * * * *

*

3

concentration: C2 ¼ 5000 mol m , flow rate: F2 ¼ 8  104 m3 s1, temperature: T2 ¼ 324 K, temperature in the jacket: Tj2 ¼ 324 K.

Finally, the customer needs a product at a fixed flow rate F ¼ 7  104 m3 s1 with the concentration C and the temperature T included between a minimum and a maximum value.

ð3Þ

*

l1 ðlj Þ is not represented inside a phase if the vector field l1 ðli ÞðxÞ ¼ x’ ¼ 0; that is if the continuous variables are constant. If m2 ðlj Þ is not represented inside a phase, it is equal to the global constraint GC: Initial and final phases are represented by double circles. If the functions set m3 is the identity (i.e. if the state variables have the same values before and after a transition occurs), it is not represented beside the arc.

In this paper, the goal is to find the optimum switching time instants which minimize a given criterion (Manon, Valentin-Roubinet, & Gilles, 2000). More details can be found in (Manon, 2001). In the following

136

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

sections, two different methods which can minimize the overall time duration of a given sequence are presented. The first synthesizes the hybrid dynamical system control by stating the optimal problem in terms of functional and constraints and the second is based on Pontryagin’s maximum principle extended to hybrid dynamical systems.

3. Optimal control design: a numerical method 3.1. Spotting the problem The process can be driven to its final phase by switching from one of its dynamic phase to another and then, allowing time to elapse until the state vector value is acceptable for the final desired hybrid region. A controllable hybrid sequence (composed of phases, each associated with its own specific time duration) which respects the constraints set, must be determined. A systematic method which can generate all the possible control sequences (Manon & Valentin-Roubinet, 1999a, b; Manon, 2001) has defined a controllable sequence, including the minimum number of locations (phases). This method is based on an analytical description of the reachable sets in the continuous state-space. It runs backward: 1. The continuous constraints in the final location define a geometrical domain in the state-space. It is analytically described. 2. Then, a possible discrete transition which reaches the final location is selected. 3. The possible values the state-vector can have before the transition occurs define a new geometrical domain called the jump-set. Its analytical description is then calculated. 4. The region of the continuous state-space which allows the system to reach the jump-set according to the dynamic of the current location is called the extended jump-set. It is analytically described. 5. The last three steps (steps 2–4) are repeated until a transition from the initial location can be fired. Then, the initial point is included into an extended jump-set. This controllability method gives a location sequence, which respects the event and the continuous constraints, i.e. a controllable sequence. But the exact switching time instants are not calculated. The transitions may be fired as soon as possible, as late as possible or somewhere in between. If they are fired ‘‘as soon as possible’’, the resulting hybrid sequence is as follows: ph ¼ fðl0 ; 0Þ; ðl1 ; 4547Þ; ðl2 ; 2228Þ; ðl1 ; 1740Þ; ðl3 ; NÞg: The times (in seconds) are the time duration the system remains in the different locations. The total time length

is then 4547+2228+1740=8515 s that is approximately 2 h 22 min. The aim is to calculate the shortest time for the system to reach a final hybrid region from an initial point included into an initial hybrid region. The state vector trajectory must respect the continuous variable constraints and the chosen discrete control sequence. Clearly, the system must switch from phase l0 to phase l1 immediately after the beginning (because x’ ¼ 0) and the system switches to the final hybrid phase l3 as soon as the continuous variables belong to m2 ðl3 Þ. Therefore, for each discrete sequence, the optimal control problem consists in calculating the intermediate switching times, i.e. in calculating the time during which the system ’ remains in each dynamic phase (where xa0). Two different methods which calculate the optimal switching times for a given sequence will be presented. But first, the existence conditions for an optimal trajectory, in the general case where the sequence is not given, are investigated. (Athans & Falb, 1966) have proven that the existence of the optimal trajectory is assumed if the number of discontinuities in the optimal trajectory is countable, i.e. if the control is piecewise continuous. Theoretically, it can happen that the state optimal trajectory runs along a constraint, which involves that the system must switch an infinite number of times between two phases (case of Zeno automaton) (Johansson, Egerstedt, Lygeros, & Sastry, 1999). But, a minimum time duration is assumed in each encountered phase for physical reasons, then the assumption that the systems which are studied are non-Zeno can be taken. With this assumption, an infinite number of transitions cannot happen in a finite time and the existence of an optimal trajectory is assumed. So, the methods presented below are applicable to any non-Zeno HDS with continuous sensors and discrete event control. In this paper, the HDS is modelled by the hybrid automaton of definition 1 with a given controllable sequence, and the existence condition of an optimal solution can also be verified. The studied class of HDS includes systems which switch between m vector fields l1 over time on the occurrence of a discrete event control S: They are called Discrete Event Controlled HDS (DECHDS). They can also be described as follows: ’ ¼ f ðxðtÞ; uðtÞÞ; xðtÞ n

ð4Þ

where xAXCR is the state vector, u ¼ ½u1 yum T AUCf0; 1gm is the control vector such that 8ta ptotb ; (iAf1; y; mg; such that ui ¼ 1 and 8jAf1; y; mg=jai; uj ¼ 0: It means that the system P is in phase i from time ta to time tb ; f ðxðtÞ; uðtÞ; tÞ ¼ m k¼1 uk :l1 ðlk Þ is a function vector of dimension n: As long as the system remains in the same phase k; the function vector f ðx; u; tÞ is equal to l1 ðlk Þ and the control uðtÞ is constant. By extension, if the controllable

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

sequence is composed of a countable number of phases, then the control uðtÞ will be piecewise constant and the system is non-Zeno when it runs from its initial value x0 to its final value xf : Thus, the existence of an optimal trajectory is verified. In this Section 3, the approach chosen to solve the optimal problem is to express it in terms of a functional and constraints. Since constraints are functions of the different vector fields l1 and of the transition functions m3 ; it is necessary to be able to calculate the vector field l1 at any time of the optimal trajectory. So, a way to achieve it will be presented. In Section 4, the PMP, usually used for continuous systems, is extended to HDSs and is applied to solve the problem. 3.2. From the HDS control to an optimal problem First, some generalities about optimization problems will be presented and then, it will be shown how the optimal control of a HDS can be designed using a classical optimization method. An optimization problem is usually modelled as follows (Ciarlet, 1990; Culioli, 1994). Given a compact subset UDVCRn and a function J : V-R; a minimum of the function J according to U must be found, i.e. an element u which verifies the following problem ðPÞ must be found: ðPÞ uAU and JðuÞ ¼ min ðJð*ÞÞ: *AU

ð5Þ

The function J is usually called functional or objective function. If the set U is equal to V; the optimization problem ðPÞ is called ‘‘without constraints’’, else, it is called ‘‘with constraints’’. A very large field of applications exists with a set U defined by 9 8 > > = < *AV= 0 : ð6Þ U¼ gi ð*Þ ¼ 0 if 1pipc ; > > ; : 0 gi ð*Þp0 if c þ 1pipc The functions gi : V-R are the constraints of the problem ðPÞ. The first (such that 1pipc0 ) defines the equality constraints and the last (such that c0 þ 1pipc) defines the inequality constraints. Now, let us see how the optimal control for a HDS can be designed. As explained in Section 3.1, the aim is to calculate the time length during which the system ’ remains in each dynamic phase (where xa0) in order to minimize the total duration of a given sequence. Let us give the definition of a dynamic sequence: Definition 2. A dynamic sequence dynseq is composed of phases with vector fields different from zero. Thus, the initial and the final phases do not appear into a dynamic sequence. The goal of the optimal problem is to minimize the overall dynamic sequence duration, therefore the

137

functional J is J¼

ns X

ð7Þ

tsi ;

i¼1

where ns is the number of phases of the dynamic sequence and tsi is the time length during which the system remains in the ith phase of the dynamic sequence, lsi Adynseq: Now, let us identify the constraints gi . The final constraints set m2 ðlf Þ is composed of cf linear inequalities. It can be directly modelled as cf inequality constraints. But, for computability reasons, it is often better to determine which final constraints may have an influence on the optimal trajectory. Therefore, the constraints m2 ðlf Þ are introduced into the optimal problem ðPÞ as two kinds: equality constraints m2eq ðlf Þ and inequality constraints m2in ðlf Þ. The problem ðPÞ is then simpler. All the other constraints m2 ðlj Þ such that lj eQ0 and lj eQf ; are constraints associated with dynamic sequence phases. They can be modelled as inequality constraints, but for the same computability reasons, it is better to determine those which have a consequence on the result (which would not be respected by an optimization without constraints) and those which have no effects. The first are included in the set m2in ðlj Þ; and the others are not modelled. The functions gi ð1pipcÞ modelling all these constraints depend on the state-vector x: Therefore, the vector fields l1 and the set of functions associated with each phase transition m3 are necessary to evaluate it. From this point, the following notations will be used: tsi is the remaining time of the control sequence in the ith phase of dynseq and xsi ¼ ½xsi;1 ; xsi;2 ; xsi;3 T is the state-vector of the ith phase of the dynamic sequence dynseq. Then, the state vector of a process in phase lk from t time unit, can be evaluated as follows: xs1 ð0Þ ¼ m3 ðl0 ; s0;s1 ; ls1 Þðx0 Þ for i ¼ 1 to ns: Z tsi xsi ðtsi Þ ¼ l1 ðlsi ÞðxðtÞÞ dt 0

xsiþ1 ð0Þ ¼ m3 ðlsi ; ssi;siþ1 ; lsiþ1 Þðxsi ðtsi ÞÞ end xk ðtÞ ¼

Z 0

t

l1 ðlk ÞðxðtÞÞ dt;

where x0 is the initial point of the trajectory. To make sure that the constraints are respected throughout the trajectory, the state-vector should be evaluated at each time instant. As a continuous evaluation of the state-vector cannot be practically realized, it is only evaluated just before and after a transition occurs. If the trajectories described by the state variables in the state-space are monotonous (cf.

138

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

Fig. 4. Monotonous vector field.

Fig. 4), the continuous evaluation does not induce any problem: indeed, if the constraints m2 ðli Þ are verified at the end of a phase, they are also verified all along this phase trajectory because the vector field remains the same. But, if the trajectory has a non-monotonous evolution (cf. Fig. 5), it is not necessarily the case. Therefore, the numerical method used to solve the optimal problem must impose small enough variations on tsi ; between two numerical iterations to be sure that, if the constraints m2 ðli Þ are violated, the state-vector at the switching time tsi is outside the constraints too. Finally, the optimal problem ðPÞ can be modelled as follows: ! ns X ðPÞ min tsi U

i¼1

with the constraints gj ðj ¼ 1; y; cÞ such that : xs1 ð0Þ ¼ m3 ðl0 ; s0;s1 ; ls1 Þðx0 Þ; for i ¼ 1 to ns  1: Z tsi xsi ðtsi Þ ¼ l1 ðlsi ÞðxðtÞÞ dt 0

xsiþ1 ð0Þ ¼ m3 ðlsi ; ssi;siþ1 ; lsiþ1 Þðxsi ðtsi ÞÞ end xsns ðtsns Þ ¼

Z

tsns 0

l1 ðlsns ÞðxðtÞÞ dt

xf ¼ m3 ðlsns ; ssns;f ; lf Þðxsns ðtsns ÞÞ and m2eq ðlf Þ ¼ 0 m2in ðlf Þp0 m2in ðlf Þp0 with i ¼ 1; y; ns : Since both the functional J and the constraints g j ð j ¼ 1; y; cÞ depend on the switching times tsi ði ¼ 1; y; ns Þ; the controllable set U changes at every new iteration of the numerical method. Therefore, the convergence of the method is not obvious and the choice of the initial guess must fit well. 3.3. Application to the product supply process An optimal control for the process modelled in Section 2 will be calculated.

Fig. 5. Non-monotonous vector field.

Two cases are considered. First, the constraints of the dynamic phases are neglected to solve the unconstrained optimization problem (UOP). Then, they are taken into account to solve the constrained optimization problem (COP). Let us present the notations used in this section and some points about the sequence to be optimized. There are an infinite number of controllable sequences to drive the process from l0 to l3 because of the loop l1  l2  l1 (see Fig. 3). The controllability method briefly described at the beginning of this Section 3 gives one of the shortest sequence (in terms of number of phases) which respects the constraints. It is equal to: seq ¼ ½l0 l1 l2 l1 l3  and the corresponding dynamic sequence is dynseq ¼ ½l1 l2 l1 . The optimal total trajectory time length will be calculated for this sequence only. Some comments on the results will prove that this leads to the absolute minimum time whatever the sequence is. 3.3.1. Preliminary remark For computability reasons, only the constraints which may have an influence on the optimal trajectory are modelled. If the evolution of the state-vector is observed from its initial point, it appears that the dynamics of the variable x2 is ‘‘slower’’ than the dynamics of x1 and x3 : Therefore, the optimal solution is achieved if the final point of the trajectory has a value for x2 as close as possible to its initial value. Since x02 ¼ 303 K and m2 ðl3 Þ includes 311px32 p315; it is deduced that x32 ¼ 311 K in the final location. 3.3.2. The unconstrained optimization problem Pu Dynseq includes three dynamic phases, thus three switching times must be calculated and the functional J is J¼

3 X

tsi ¼ ts1 þ ts2 þ ts3 :

ð8Þ

i¼1

The final value of the state-vector is x3 ¼ ½x31 ; 311; x33 T and the constraints for the final phase are m2 ðl3 Þ: 2200px31 p2300; 311px32 p315; 0px33 p10. There-

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

Continuous state vector in phase ls3 ¼ l1 at time t ¼ ts3

fore, the equality constraint is x3;2 ð0Þ ¼ 311:

ð9Þ

2200px3;1 ð0Þp2300; 0px3;3 ð0Þp10;

¼ 0s :

2

2300

3

6 7 xs3 ðts3 Þ ¼ 4 311 5: 9:07

The four inequality constraints are ð10Þ

The resulting optimal problem Pu can now be explicitly given as Pu : minðts1 þ ts2 þ ts3 Þ g1 : x3;2 ð0Þ  311 ¼ 0; g2 : x3;1 ð0Þ  2300p0; g3 : x3;1 ð0Þ þ 2200p0; g4 : x3;3 ð0Þ  10p0; g5 : x3;3 ð0Þp0; with x0 ¼ ½2500 303 10T ; xs1 ð0Þ ¼ x0 and xs1 ðts1 Þ ¼

Z

ts1

l1 ðl1 ÞðxðtÞÞ dt; Z ts2 l1 ðl2 ÞðxðtÞÞ dt; xs2 ð0Þ ¼ xs1 ðts1 Þ and xs2 ðts2 Þ ¼ 0 Z ts3 l1 ðl3 ÞðxðtÞÞ dt: xs3 ð0Þ ¼ xs2 ðts2 Þ and xs3 ðts3 Þ ¼ 0

The optimal state trajectory solution of the UOP is presented in Fig. 6 as well as the solution in COP, which is explained later. The evolution of the state-variables over time are also given in Figs. 7–9 in both cases. It appears that the third transition, which is compulsory to get a controllable sequence so that the GCs are respected, is unnecessary in UOP (ts3 ¼ 0). It means that the same result would be obtained if the sequence to optimize in the UOP was sequ ¼ ½l0 l1 l2 l3  and the corresponding dynamic sequence dynsequ ¼ ½l1 l2 : The minimum global time length is then ts1 þ ts2 þ ts3 ¼ 6986 s ¼ 1 h 56 min. This is hopefully less than with the ‘‘as soon as’’ transition strategy (2 h 22 min) explained at the beginning of this Section 3. 3.3.3. The constrained optimization problem, Pc The dynamic sequence to optimize is still the same, dynseqc ¼ ½l1 l2 l1 : There are six inequality constraints

0

It is solved by a program based on the instruction ‘‘constr’’ of the Matlab optimization toolbox. This instruction finds the constrained minimum of a function of several variables. The functions g1 to g5 are calculated at the transition instants owing to the different vector fields l1 and the different transition functions m3 : The result is Continuous state vector in phase l0 at time t ¼ 0 : 2 3 2500 6 7 xð0Þ ¼ 4 303 5:

Fig. 6. The optimal state trajectory.

10 Continuous state vector in phase ls1 ¼ l1 at time t ¼ ts1 ¼ 5434 s :

2

1900

3

6 7 xs1 ðts1 Þ ¼ 4 306 5: 8:91 Continuous state vector in phase ls2 ¼ l2 at time t ¼ ts2 ¼ 1552 s : 2 3 2300 6 7 xs2 ðts2 Þ ¼ 4 311 5: 9:07

139

Fig. 7. Concentration versus time.

140

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

Fig. 8. Temperature versus time.

Fig. 9. Volume versus time.

per phase (a minimum and a maximum value for each state-variable) and the dynamic sequence is composed of three phases, so eighteen new inequality constraints could be added to the unconstrained optimal problem Pu to get the constrained optimal problem Pc . But, the unconstrained optimal problem solution highlights that only the constraint on the concentration at the end of phase ls1 ¼ l1 can be violated with the sequence dynseq. Therefore, only one inequality constraint is added.

The constrained optimal problem Pc is solved with the same numerical tool as previously. The result is

xs1;1 X2000:

ð11Þ

Then, the constrained optimal problem Pc is given below: Pc : minðts1 þ ts2 þ ts3 Þ

Continuous state vector in phase ls1 ¼ l1 at time t ¼ ts1 ¼ 4460 s :

2

3 2000 6 7 xs1 ðts1 Þ ¼ 4 305:6 5: 9:11

Continuous state vector in phase ls2 ¼ l2 at time t ¼ ts2 2486

3

6 7 xs2 ðts2 Þ ¼ 4 312 5: 9:31

g2 : x3;1 ð0Þ  2300p0; g3 : x3;1 ð0Þ þ 2200p0;

Continuous state vector in phase ls3 ¼ l1 at time t ¼ ts3 ¼ 1484 s : 2 3 2300 6 7 xs3 ðts3 Þ ¼ 4 311 5:

g4 : x3;3 ð0Þ  10p0; g5 : x3;3 ð0Þp0; g6 : xs1;1 ð0Þ þ 2000p0;

9:01

with

xs1 ð0Þ ¼ x0 and xs1 ðts1 Þ ¼

10

¼ 2033 s : 2

g1 : x3;2 ð0Þ  311 ¼ 0;

x0 ¼ ½2500 303 10T ;

Continuous state vector in phase l0 at time t ¼ 0 : 2 3 2500 6 7 xð0Þ ¼ 4 303 5:

Z

ts1

l1 ðl1 ÞðxðtÞÞ dt; Z ts2 l1 ðl2 ÞðxðtÞÞ dt; xs2 ð0Þ ¼ xs1 ðts1 Þ and xs2 ðts2 Þ ¼ 0 Z ts3 l1 ðl3 ÞðxðtÞÞ dt: xs3 ð0Þ ¼ xs2 ðts2 Þ and xs3 ðts3 Þ ¼ 0

0

The optimal trajectory and the evolution of the statevariables are presented in Figs. 6–9. In the constrained optimization problem, all the dynamic phases of the sequence are necessary because no time length tsi is null (hopefully as stated by the controllability analysis). The total duration time of the sequence seqc ¼ ½l0 l1 l2 l1 l3  is then 4460+2033+1484=7977 s, that is approximately 2 h 13 min. This is faster than if the transitions are fired as soon as possible (2 h 22 min) and longer than in the UOP solution. This seems logical since the more the

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

system is constrained the more time it may need to reach the final phase.

4.1. The PMP extended to HDS The PMP can be applied to any system with a piecewise continuous control (Pontryagin, Boltyanskii, Gamkrelidze, & Mishchenko, 1964). It was proven at the beginning of Section 3 that this condition is true for any non-Zeno automaton, in particular for the product supply process model. Let us give some recalls about Pontryagin’s theory. First, let us define a functional J; used as a performance criterion: J¼

tf

f0 ðx; uÞ dt;

ð12Þ

t0

where f0 ðx; uÞ is a given function. The optimal problem can be stated as follows: Let us consider two points x0 AX and x f AX; the problem is to find, if it exists, the control uðtÞAU which drives the system from x0 to x f and which minimizes the functional J: If the function f0 ðx; uÞ ¼ 1; the functional J becomes J ¼ tf  t0 and the control uðtÞ is optimal if the tim e which is necessary to drive the system from x0 to x f is minimum. This is the problem which will be dealt in this paper and it is called the time optimal control problem. Let us define a new state variable x0 to translate the optimal control problem into a system of algebrodifferential equations. In fact, because of the criterion f0 ðx; uÞ ¼ 1; x0 represents the time. dx0 ¼ f0 ðx; uÞ: dt

ð14Þ

where

In this section, it will be shown how an optimal control problem for HDS with a discrete event control (DECHDS) can be translated into a system of algebraic and differential equations. This translation will be achieved owing to PMP which gives necessary conditions for a control to be optimal. First, some recalls about the PMP will be given. Next, a hybrid version of the PMP which takes into account the controlled and autonomous transitions will be given. The notations are the same as for hybrid automaton, so it will be directly possible to apply this theory to the product supply process described in Section 2. This will be achieved in Section 4.2.

Z

The differential equation system (4) can be augmented as follows by adding dx0 =dt: dx* * uÞ; ¼ fðx; dt

4. Optimal control with the Pontryagin’s maximum principle

141

ð13Þ

nþ1 * x* ¼ ðx0 ; xÞAXCR and " # f0 ðx; uÞ ¼ 1 nþ1 * uÞ ¼ * fðx; : AXCR f ðx; uÞ

Then, the optimal control problem can be stated as follows: * Let us consider two points ð0; x0 ÞAX* and ðx0f ; x f ÞAX: It must be found, if it exists, the control uðtÞAU which drives the system from ð0; x0 Þ to ðx0f ; x f Þ such that x0f is minimum (cf. Fig. 10). Now, let us define the adjoint vector (Lagrange multipliers): pAPCRn and the augmented adjoint vector nþ1 * : p* ¼ ðp0 ; pÞAPCR The following differential equations system can be introduced: * uÞT dp* qfðx; ¼ p: * ð15Þ dt qx* Systems (14) and (15) can be linked by the function H: * uÞ: Hðx; p; * uÞ ¼ p*T fðx; ð16Þ It is easily verified that systems (14) and (15) can be written in the Hamiltonian system form: dx* qH ¼ ; ð17Þ dt qp* dp* qH ¼ : ð18Þ dt qx* Thus, H is the Hamiltonian function for the increased systems (14) and (15). In case of time optimal control problem, Eq. (16) can be written as follows: Hðp0 ; p; x; uÞ ¼ p0 þ pT f ðx; uÞ:

Fig. 10. The optimal control.

ð19Þ

142

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

4. For a given set ðp0 ; p; x; uÞ at time t; if t is not a switching time and if the final instant tf is not fixed:

For hybrid models, this equation is Hðp0 ; p; x; uÞ ¼

m X

Hi ðp0 ; p; xÞui ðtÞ;

ð20Þ

8tA½t0 ; tf ; Hðp0 ; p; x; uÞ ¼ 0:

ð24Þ

i¼1

where Hi is the Hamiltonian function associated with phase li . In the hybrid automaton model, f is defined by the dynamic l1 ðli Þ: Then, Eq. (20) can be written as below Hi ðp0 ; p; xÞ ¼ p0 þ

n X

pj l1 ðli Þðxj Þ:

ð21Þ

j¼1

Finally, conditions about the extremity of the trajectory must be added to describe completely the problem. When the state-vector optimal trajectory reaches its final value, it is considered that kf inequality constraints among cf defining m2 ðlf Þ are equal to zero. These kf equalities describe a manifold of dimension n  kf and then 1pkf pminðn; cf Þ. Let us now define the plan Tf which is tangent to the previous manifold at point x f : If the control uðtÞ and the corresponding trajectory xðtÞ are optimal, then the adjoint vector p satisfies the transversality conditions at the final extremity of the trajectory xðtÞ; i.e. pðtf Þ is orthogonal to Tf : Formally,   qT f T pðtf Þ ¼ pf jt¼tf ; ð22Þ qx where pf is a kf length vector. PMP adapted to the continuous parts of the DECHDS can now be presented. It gives necessary conditions for a control to be optimal. A new version of the PMP extended to both parts (continuous and discrete) of the DECHDS will be presented later.

5. At the final time tf ; the vector p must verify the transversality conditions stated by Eq. (22). It is noticed that, when a transition from one phase to another occurs, the adjoint vector may not be continuous. It is continuous everywhere else, so it is said to be ‘‘piecewise continuous’’. For the same reason, condition one is true ‘‘almost everywhere’’ because it may be false when a transition occurs. Now, let us see how to extend PMP to Discrete Event Controlled HDS (DECHDS). When the DECHDS switches from dynamics l1 ðli Þ to dynamics l1 ðlj Þ at time tsi ; the discrete transition trij can be either controlled or autonomous. If one of the continuous state variables overpasses a limit of the constraint set m2 ðli Þ after t ¼ tsi ; the system must leave the current phase li at t ¼ tsi ; to another phase lj in order to respect the constraints. Therefore, the switch is imposed by the plant. sij is an uncontrolled event and the transition is autonomous. If the continuous variables are not on a limit of the constraint set m2 ðli Þ at t ¼ tsi ; sij is a controlled event and the transition is also controlled. If a controlled transition from phase li to phase lj occurs in a DECHDS, new equations are issued of Eq. (23). The Hamiltonian function Hi at phase li and the Hamiltonian function Hj at phase lj respect the following conditions: sup H ¼ Hi > Hj at t ¼ tsi  dt and uAU

Hi oHj ¼ sup H at t ¼ tsi þ dt: uAU

Then sup H ¼ Hi ¼ Hj at t ¼ tsi :

ð25Þ

uAU

4.1.1. The Pontryagin’s maximum principle (continuous parts) If uAU is an admissible (piecewise continuous) optimal control and xAX the associated optimal trajectory for problems (4) and (12), then the scalar value p0 and the piecewise continuous vector function pAP are defined such that: 1. The set ðp0 ; p; x;uÞ satisfies the Hamiltonian system almost everywhere. 2. The scalar value p0 is negative constant and can be equal to 1. 3. For a given set ðp0 ; p; x; uÞ at time t; the following maximal condition is verified if t is not a switching time from a phase li to a phase lj : Hðp0; p; x; uÞ ¼ sup

m X

uAU i¼1

Hi ðp0 ; p; xÞui ðtÞ:

ð23Þ

Then transversality conditions imply: pðtsi þ dtÞ ¼ pðtsi  dtÞ:

ð26Þ

If an autonomous transition from phase li to phase lj occurs in a DECHDS, ki inequality constraints among ci which define m2 ðli Þ are equal to zero. These ki equalities define a manifold of dimension n  ki . Then, the plan T i which is tangent to this manifold at the point xi reached by the state-vector can be defined and the next equations can be written as   qT i T Hj ¼ Hi þ pijt¼tsi : ð27Þ qt  pðtsi þ dtÞ ¼ pðtsi  dtÞ 

qT i qx

T pijt¼tsi ;

ð28Þ

where pi is a ki length vector. Eqs. (25)–(28) are transversality conditions. Mathematically rigorous demonstration for Eqs. (25)–(28) can be found in (Bryson

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

& Ho, 1969) in terms of continuous systems with statevector constraints. (Riedinger, Zanne, & Kratz, 1999) and (Riedinger, 2000) have also given a formal demonstration for hybrid systems. Moreover, they have adapted these equations to systems which involve a state-vector discontinuity when a discrete transition occurs. xðtsi þ dtÞ ¼ m3 ðli ; sij ; lj Þðxðtsi  dtÞÞaxðtsi  dtÞ

ð29Þ

Eqs. (25) and (27) lead to the following equations: *

For a controlled transition: Hi ðp0 ; p; xÞ ¼ Hj ðp0 ; p; m3 ðli ; sij ; lj ÞðxÞÞ:

*

ð30Þ

For an autonomous transition: Hj ðp0 ; p; m3 ðli ; sij ; lj ÞðxÞÞ   qT i T ¼ Hi ðp0 ; p; xÞ þ pijt¼tsi : qt

ð31Þ

Finally, a hybrid version of the PMP can be stated: 4.1.2. The Pontryagin’s maximum principle extended to DECHDS The necessary conditions 1–5 of the PMP adapted to the continuous parts of the DECHDS must be respected and the sixth one must be verified: at the switching time from phase li to phase lj , t ¼ tsi , transversality conditions are verified: *

*

If the transition is controlled, Eqs. (25) (or (30)) and (26) are true, if the transition is autonomous, Eqs. (27) (or (31)) and (28) are true.

The fourth condition of the PMP states that the Hamiltonian function H in a given phase li is constant. This value is chosen equal to zero in the first dynamic phase ls1 : If the transition to the second phase ls2 is controlled, then Eq. (25) involves that the Hamiltonian function remains equal to zero. But if this transition is autonomous, then the Hamiltonian function is shifted by a value equal to ½qT i =qtT pijt¼tsi (cf. Eq. (27)). The same variation can be observed at any transition of the control sequence. This shifting means physically that the process would be driven faster to its final phase if the transition would occur later but the transition happens in order to respect the constraints. As with the numerical method, two cases are considered. First, the constraints of the dynamic phases are neglected to solve the UOP. Then, they are taken into account to solve the COP. With this method, the solution of the UOP will be necessary to calculate the optimal switching times of the COP. Indeed, in the UOP, every transition is controlled, thus, Eqs. (25) and (26) systematically belong to the system to solve. If the corresponding optimal trajectory remains in the global state constraints GC; then the optimal trajectory in

143

the COP is the same. If not, Eqs. (27) and (28) for the transitions not included into GC; must be added to the UOP representation (instead of Eqs. (25) and (26)). The subject of this section was to explain how an optimal control problem for DECHDS can be translated into a system of algebraic and differential equations. This translation is achieved thanks to the PMP theory which is extended to hybrid systems. The PMP gives necessary but no sufficient conditions. This points out the problems of existence and uniqueness of an optimal control. By now, it has only been proven that an optimal control exists and is unique for linear systems (Pontryagin, 64). As the system of equations to be solved in this paper can be non-linear, the solution may not be unique and computational problems may also appear. 4.2. Application to the product supply process The time optimal control for the process modelled in Section 2 will be calculated with this second method. First, some points about the sequence to be optimized and about a possible final point for the state variables are discussed. 4.2.1. Some preliminaries In Section 3.3, it was explained that the absolute minimum time in COP is obtained by optimizing the shortest controllable sequence: seqc ¼ ½l0 l1 l2 l1 l3 : The corresponding dynamic sequence is dynseqc ¼ ½l1 l2 l1 : Since the PMP optimization method cannot model inequality constraints, a possible final point x3 ; must be defined. The final point x3 has three components: x3 ¼ ½x31 x32 x33 T : In Section 3.3, it was established that x32 ¼ 311 K. Some analysis could lead to the possible values for x31 and x33 to get the optimal solution. But, in order to avoid this analysis, the results of the Section 3 are used. Therefore, x31 ¼ 2300 mol m3 and x33 is free. So, the preliminary conclusions are: *

The final point x3 must be equal to: x3 ¼ ½2300 311 x33 T with 0px33 p10:

*

ð32Þ

The sequence to be optimized and the corresponding dynamic sequences are: seq ¼ ½l0 l1 l2 l1 l3  dynseq ¼ ½l1 l2 l1 :

ð33Þ

4.2.2. The unconstrained optimization problem The system of algebraic and non-linear differential equations in the UOP is made of 24 equations established as follows: Phase ls1 ¼ l1 : At time t ¼ 0; the process is in dynamic phase l1 : Eq. (21) and condition 2 of the PMP lead to: H1 ¼ 1 þ pTs1 l1 ðl1 ÞðxÞ.

144

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

From Eq. (24), H1 ¼ 0 whatever tpts1 and at t ¼ 0; x0 is fixed, so the first equation (Eq. (A.1)) is obtained: H1 ¼ 1 þ pTs1 l1 ðl1 Þ ¼ 0 at t ¼ 0: During phase l1 ; the dynamic is l1 ðl1 Þ. Therefore, three differential equations (Eqs. (A.2)–(A.4)) must be added: x’ s1 ¼ l1 ðl1 Þ: From Eq. (18), three differential equations (Eqs. (A.5)–(A.7)) identify the adjoint vector: dp*s1 =dt ¼ ðqH1 =qxÞ * T: Switch from phase ls1 ¼ l1 to phase ls2 ¼ l2 : At time t ¼ ts1 ; the process switches from l1 to l2 : In the UOP, all events are controlled because they can happen at any time, without taking care of any constraints. Thus, function m3 ðl1 ; s12 ; l2 ÞðxÞ and Eq. (26) respectively assume the continuity of the state-vector and of the adjoint vector when the process switches xs2 ð0Þ ¼ xs1 ðts1 Þ and ps2 ð0Þ ¼ ps1 ðts1 Þ: Eq. (21) and condition 2 of the PMP lead to: H2 ¼ 1 þ pTs2 l1 ðl2 ÞðxÞ: From Eq. (25), a new equation is found (Eq. (A.8)): H2  H1 ¼ 0: Phase ls2 ¼ l2 : the state and the adjoint vector calculation in phase l2 gives six new differential equations (Eqs. (A.9)–(A.14)) x’ s2 ¼ l1 ðl2 Þ and dp*s2 =dt ¼ ðqH2 =qxÞ * T: Switch from phase ls2 ¼ l2 to phase ls3 ¼ l1 : At time t ¼ ts2 ; the process switches from l2 to l1 : As for the first switch, it is a controlled event. Function m3 ðl2 ; s21 ; l1 ÞðxÞ and Eq. (A.26) respectively, assume the continuity of the state-vector and of the adjoint vector when the process switches: xs3 ð0Þ ¼ xs2 ðts2 Þ and ps3 ð0Þ ¼ ps2 ðts2 Þ: Eq. (21) and condition 2 of the PMP lead to: H1 ¼ 1 þ pTs3 l1 ðl1 ÞðxÞ: From Eq. (A.25), a new equation is found (Eq. (A.15)): H1  H2 ¼ 0: Phase ls3 ¼ l1 : the state and the adjoint vector calculation in phase l1 gives six new differential equations (Eqs. (A.16)–(A.21)) x’ s3 ¼ l1 ðl1 Þ and dp*s3 =dt ¼ ðqH1 =qxÞ * T: At the final time t ¼ tf ¼ ts1 þ ts2 þ ts3 : the continuous state variables are such that k3 inequality constraints of the final constraint region m2 ðl3 Þ are equal to zero. According to the preliminary conclusions, the optimal trajectory ends on the following straight-line x3 ¼ ð2300; 311; x33 Þ: Therefore, k3 ¼ 2 and two new equations are defined (Eqs. (A.22) and (A.23)). From Eq. (22), it is deduced that ps2 ðts2 Þ ¼ ½p31 p32 0T : Only the third equation (Eq. (A.24)) is interesting. In total, 24 equations, explicitly given in Appendix A.1, are obtained. They are algebraic or differential and some are non-linear. There are 24 unknown variables:

vector ps1 ð0Þ; time ts1 ; vectors ps1 ðts1 Þ and xs1 ðts1 Þ; time ts2 ; vectors ps2 ðts2 Þ and xs2 ðts2 Þ; time ts3 ; vectors ps3 ðts3 Þ and xs3 ðts3 Þ: At each iteration, function ‘‘ode45’’ solves the differential equations and then function ‘‘fsolve’’ of the matlab optimization toolbox is used to solve numerically the system of equations. The result is Phase l0 at time t ¼ 0: 2 3 2 3 2500 5:14 6 7 6 7 x0 ð0Þ ¼ 4 303 5 p0 ð0Þ ¼ 4 308 5: 10

699

Phase ls1 ¼ l1 at time t ¼ ts1 ¼ 5434 s: 2 3 2 3 1900 6:85 6 7 6 7 xs1 ðts1 Þ ¼ 4 306 5 ps1 ðts1 Þ ¼ 4 766 5: 8:91 174 Phase ls2 ¼ l2 at time t ¼ ts2 ¼ 1552 s: 2 3 2 3 2300 7:86 6 6 7 7 xs2 ðts2 Þ ¼ 4 311 5 ps2 ðts2 Þ ¼ 4 1062 5: 9:07

0

Phase ls3 ¼ l1 at time t ¼ ts3 ¼ 0 s: 2 3 2 3 2300 7:86 6 7 6 7 xs3 ðts3 Þ ¼ 4 311 5 ps3 ðts3 Þ ¼ 4 1062 5: 9:07

0

It can be verified that the switching times (ts1 ; ts2 and ts3 ) and the corresponding state-vectors have the same values as those calculated with the first method (cf. Section 3.3.1). The optimal state trajectories are already presented in Figs. 6–9. The optimal adjoint trajectory and the evolution of the Hamiltonian function are given on Figs. 11 and 12.

Fig. 11. The optimal adjoint trajectory.

145

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

and xs2 ðts2 Þ; time ts3 ; vectors ps3 ðts3 Þ and xs3 ðts3 Þ. The solution is calculated with the same algorithm than in the UOP. The resulting values are: Phase l0 at time t ¼ 0 s: 2 3 2 3 2500 4:63 6 7 6 7 xð0Þ ¼ 4 303 5 pð0Þ ¼ 4 391 5: 10

668

Phase ls1 ¼ l1 at time t ¼ ts1 ¼ 4460 s: 2 3 2 3 2000 5:85 6 7 6 7 xs1 ðts1 Þ ¼ 4 305:6 5 ps1 ðts1 Þ ¼ 4 822 5p11 ¼ 3:54: 9:11 244 Fig. 12. Hamiltonian function over time.

In the UOP, as all the transitions are controlled, it can be verified that the Hamiltonian function is equal to zero all along the trajectory (cf. Fig. 12).

4.2.3. The constrained optimization problem The dynamic sequence to optimize is still dynseqc ¼ ½l1 l2 l1 . It must be determined if the events s12 and s21 which lead to dynseqc are controlled or not to deduce the resulting equations of the system. In UOP, the optimal trajectory of the state variables verifies the dynamic phases constraints everywhere except around the first transition from l1 to l2 ; where the concentration is below the allowed minimum: 2000 mol m3 (cf. Figs. 6 and 7). Therefore, in the constrained optimization problem, the process must switch from l1 ðl1 Þ to l1 ðl2 Þ when x1 ¼ 2000: It defines an autonomous transition on the occurrence of s12 . In the first guess, the event s21 will be considered as controlled. If the result does not respect the constraints, a new calculation will be achieved with s21 as an uncontrolled event. The system of non-linear algebraic and differential equations is established by the same method as in the UOP. The equations due to the autonomous transition replace Eq. (A.8) of the previous system of equations given in Appendix A.1. Then, at time t ¼ ts1 ; x1 ¼ 2000: It defines the plan T 1 and gives the new equation Eq. (A.32) in Appendix A.2. Eq. (28) gives the next three equations: ps2 ð0Þ ¼ ps1 ðt1 Þ  ½p11 0 0T ; and leads to Eq. (A.33) in Appendix A.2. Eq. (27) gives a new equation (Eq. (A.34) in Appendix A.2): H2 ¼ H1 þ ðqT 1 =qtÞT p1jt¼t : Then, the overall system is composed of 26 equations explicitly given in Appendix A.2. There are 26 unknown variables: vector ps1 ð0Þ; time ts1 ; vectors ps1 ðts1 Þ and xs1 ðts1 Þ; shifting value p11 on the adjoint vector at time t ¼ ts1 and initial value ps2;1 ð0Þ; time ts2 ; vectors ps2 ðts2 Þ

Phase ls2 ¼ l2 at time t ¼ ts2 ¼ 2033 s: 2 3 2 3 2486 11:2 6 7 6 7 xs2 ðts2 Þ ¼ 4 312 5 ps2 ðts2 Þ ¼ 4 1247 5: 9:31

101

Phase ls3 ¼ l1 at time t ¼ ts3 ¼ 1484 s: 2 3 2 3 2300 12:2 6 7 6 7 xs3 ðts3 Þ ¼ 4 311 5 ps3 ðts3 Þ ¼ 4 1612 5: 9:01 0 Once again, it can be verified that the switching times (ts1 ; ts2 and ts3 ) and the corresponding state-vectors have the same values as those calculated with the numerical method (cf. Section 3.3.2). The optimal trajectories are presented in Figs. 6–9, 11 and 12. The shifting value p11 is noted as pi in Fig. 11. It can be noticed in Fig. 12 that, in the COP, the Hamiltonian function is discontinuous at t ¼ ts1 because an autonomous transition occurs. The adjoint trajectory is also discontinuous at that time (cf. Fig. 11).

5. Conclusion and perspectives In this paper, two different methods which calculate the time optimal control for DECHDS, i.e. HDSs with a piecewise continuous control were presented. In Section 2, the hybrid automaton model was presented and applied to a product supply process. In Section 3, the first optimization method which translates the time optimal control problem for DECHDS in terms of functional and constraints was explained. First, the problem position was clarified and the existence conditions of an optimal solution was investigated. Then, the translation in terms of functional and constraints was presented. Finally, this method was applied to the product supply process. In Section 4, the second optimal control synthesis method which extends the PMP to HDS was developed. First, it was verified that the PMP could be applied to HDS with a discrete event control. Then, the necessary

146

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

conditions to transform the optimal control problem into an algebraic and differential equations system was presented. Then, as in the previous section, the method was applied to the same process. Hopefully, the results were the same. Both methods can deal with HDS with non-linear vector fields. They have no limitations with respect to: *

the number of phases or discrete states: Q ¼ flj ; jAf1; y; mgg with mAN;

*

The authors thank Franc-oise Couenne, associate researcher at the French Scientific Researches National Center (CNRS), for her remarks and comments about the topic. They also thank the National club EEA/SEE work group on HDS and the GDR Automatique for valuable discussions on these topics.

the number of continuous variables: xAVDRn with nAN;

*

Acknowledgements

the length of the sequence to optimize: jseqjAN:

The second method, which is based on a deep mathematical theory (the PMP), is nevertheless bulkier to use practically. With this method, the two transitions types (autonomous and controlled) do not lead to the same equations and the type of each transition is not initially known. Therefore, a first choice must be made without going into details. With the first method, this problem is irrelevant. Therefore, a possible design method may be to find the optimal control with the first method and to validate it with the second one. The optimization methods that have been presented here calculate the optimal control for a given sequence. It must be remembered that (Bonnard & Pelletier, 1995) have demonstrated that the optimal control of a nonlinear continuous system of dimension 2 in an UOP accepts a bang–bang optimal control with at most one switching. The same kind of studies could be led to determine which sequence to choose in the constrained optimization problem. The synthesis method does not impose a maximum number of continuous state variables. Nevertheless, the number of equations to solve increases rapidly, so only a numerical solution is calculated. But the possible nonlinearity of the dynamics introduces both convergence and uniqueness of solution problems (local vs. global minimum). Therefore, the initial guess for the numerical solution calculation must be adequately chosen. A possible method to identify an acceptable guess is to linearize the system around a point and to calculate the optimal solution which is then defined and unique. This first approximation can be used as the initial guess for the calculation with the non-linear model. In addition to these numerical convergence problems, it is recalled that the existence and the uniqueness of an optimal solution is not systematic, as explained at the end of Section 4.1. These last two points are classical problems of nonlinear control and more works to clarify it for DECHDS should be interesting to realize.

Appendix A A.1. The equations system in the unconstrained optimization problem with the PMP At t ¼ 0 in phase ls1 ¼ l1 ; x ¼ xs1 ð0Þ ¼ x0 ¼ ½x0;1 x0;2 x0;3 T    F1 ðC1  x0;1 Þ H1 ð0Þ ¼ 1 þ p1;1 x0;3   F1 ðT1  x0;2 Þ UAWm ðTj1  x0;2 Þ þ p1;2 þ x0;3 dCp x0;3  ¼ 0: ðA:1Þ þp1;3 ðF1  FÞ jt¼0

During phase ls1 ¼ l1 ; 0ptpts1 ; x ¼ xs1 ¼ ½xs1;1 xs1;2 xs1;3 T and p ¼ ps1 ¼ ½ps1;1 ps1;2 ps1;3 T : x’ s1;1 ¼

F1 ðC1  xs1;1 Þ ; xs1;3

ðA:2Þ

x’ s1;2 ¼

F1 ðT1  xs1;2 Þ UAWm ðTj1  xs1;2 Þ þ ; xs1;3 dCp xs1;3

ðA:3Þ

x’ s1;3 ¼ F1  F; p’s1;1 ¼

p’s1;2 ¼

p’s1;3 ¼

F1 ps1;1 ; xs1;3  F1 UAWm þ ps1;2 ; xs1;3 dCp xs1;3

ðA:4Þ

ðA:5Þ



ðA:6Þ

F1 ðC1  xs1;1 Þ ps1;1 x2s1;3 " # F1 ðT1  xs1;2 Þ UAWm ðTj1  xs1;2 Þ þ þ ps1;2 ; x2s1;3 dCp x2s1;3 ðA:7Þ

At the controlled transition from ls1 to ls2 ; x ¼ xs1 ðts1 Þ ¼ xs2 ð0Þ ¼ ½xts1;1 xts1;2 xts1;3 T and p ¼ ps1 ðts1 Þ¼ ps2 ð0Þ ¼ ½pts1;1 pts1;2 pts1;3 T

147

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

Hs2 ð0Þ  Hs1 ðts1 Þ   F2 ðC2  xts1;1 Þ  F1 ðC1  xts1;1 Þ ¼ pts1;1 xts1;3  F2 ðT2  xts1;2 Þ  F1 ðT1  xts1;2 Þ þ pts1;2 xts1;3  UAWm ðTj2  Tj1 Þ þ dCp xts1;3 þ pts1;3 ðF2  F1 Þ ¼ 0:

F1 ps3;1 xs3;3   F1 UAWm ps3;2 ¼ þ xs3;3 dCp xs3;3

p’s3;1 ¼

ðA:19Þ

p’s3;2

ðA:20Þ

p’s3;3 ¼ ðA:8Þ

F1 ðC1  xs3;1 Þ ps3;1 x2s3;3 " # F1 ðT1  xs3;2 Þ UAWm ðTj1  xs3;2 Þ þ þ ps3;2 : x2s3;3 dCp x2s3;3

During phase ls2 ¼ l2 ; 0ptpts2 ; x ¼ xs2 ¼ ½xs2;1 xs2;2 xs2;3 T and p ¼ ps2 ¼ ½ps2;1 ps2;2 ps2;3 T : x’ s2;1 ¼ x’ s2;2

F2 ðC2  xs2;1 Þ xs2;3

F2 ðT2  xs2;2 Þ UAWm ðTj2  xs2;2 Þ ¼ þ xs2;3 dCp xs2;3

x’ s2;3 ¼ F2  F F2 ps2;1 xs2;3   F2 UAWm ¼ þ ps2;2 xs2;3 dCp xs2;3

ðA:9Þ ðA:10Þ ðA:11Þ

p’s2;1 ¼

ðA:12Þ

p’s2;2

ðA:13Þ

p’s2;3 ¼

F2 ðC2  xs2;1 Þ ps2;1 x2s2;3 " # F2 ðT2  xs2;2 Þ UAWm ðTj2  xs2;2 Þ þ þ ps2;2 : x2s2;3 dCp x2s2;3 ðA:14Þ

At the controlled transition from ls2 to ls3 ; x ¼ xs2 ðts2 Þ¼ xs3 ð0Þ ¼ ½xts2;1 xts2;2 xts2;3 T and p ¼ ps2 ðts2 Þ¼ ps3 ð0Þ ¼ ½pts2;1 pts2;2 pts2;3 T Hs3 ð0Þ  Hs2 ðts1 Þ   F1 ðC1  xts2;1 Þ  F2 ðC2  xts2;1 Þ ¼ pts2;1 xts2;3  F1 ðT1  xts2;2 Þ  F2 ðT2  xts2;2 Þ þ pts2;2 xts2;3  UAWm ðTj1  Tj2 Þ þ dCp xts2;3 þ pts2;3 ðF1  F2 Þ ¼ 0:

ðA:15Þ

0ptpts3 ; x ¼ xs3 ¼ During phase ls3 ¼ l1 ; ½xs3;1 xs3;2 xs3;3 T and p ¼ ps3 ¼ ½ps3;1 ps3;2 ps3;3 T : x’ s3;1

F1 ðC1  xs3;1 Þ ¼ xs3;3

x’ s3;2

F1 ðT1  xs3;2 Þ UAWm ðTj1  xs3;2 Þ ¼ þ xs3;3 dCp xs3;3

x’ s3;3 ¼ F1  F

ðA:16Þ ðA:17Þ ðA:18Þ

ðA:21Þ in phase ls3 ¼ l1 ; x ¼ xs3 ðts3 Þ ¼ At t ¼ ts3 ½xs3;1 xs3;2 xs3;3 T and p ¼ ps3 ðts3 Þ ¼ ½ps3;1 ps3;2 ps3;3 T xs3;1 ¼ 2300:

ðA:22Þ

xs3;2 ¼ 311: 2 qðxs3;1  2300Þ 6 qxs3;1 6 6 qðx  2300Þ s3;1 6 ps3 ¼ 6 6 qxs3;2 6 4 qðxs3;1  2300Þ qxs3;3 2 3 ps3;1 6 7 ¼ 4 ps3;2 5:

ðA:23Þ 3 qðxs3;2  311Þ 7 qxs3;1 7" # qðxs3;2  311Þ 7 7 ps3;1 7 7 ps3;2 qxs3;2 7 qðxs3;2  311Þ 5 qxs3;3

0 gives the following equation: Ps3;3 ¼ 0:

ðA:24Þ

A.2. The equations system in the constrained optimization problem with the PMP At t ¼ 0 in phase ls1 ¼ l1 ; x ¼ xs1 ð0Þ ¼ x0 ¼ ½x0;1 x0;2 x0;3 T    F1 ðC1  x0;1 Þ H1 ð0Þ ¼ 1 þ p1;1 þ p1;2 x0;3   F1 ðT1  x0;2 Þ UAWm ðTj1  x0;2 Þ þ x0;3 dCp x0;3  þp1;3 ðF1  FÞ ¼ 0: ðA:25Þ jt¼0

0ptpts1 ; x ¼ xs1 ¼ During phase ls1 ¼ l1 ; ½xs1;1 xs1;2 xs1;3 T and p ¼ ps1 ¼ ½ps1;1 ps1;2 ps1;3 T : x’ s1;1 ¼

F1 ðC1  xs1;1 Þ xs1;3

ðA:26Þ

x’ s1;2 ¼

F1 ðT1  xs1;2 Þ UAWm ðTj1  xs1;2 Þ þ xs1;3 dCp xs1;3

ðA:27Þ

148

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

x’ s1;3 ¼ F1  F

ðA:28Þ

F1 ps1;1 xs1;3   F1 UAWm ¼ þ ps1;2 xs1;3 dCp xs1;3

p’s1;1 ¼ p’s1;2 p’s1;3

ðA:29Þ p’s2;2 ¼

F2 ps2;1 xs2;3  F2 UAWm þ ps2;2 xs2;3 dCp xs2;3

ðA:38Þ



ðA:39Þ

ðA:30Þ

F1 ðC1  xs1;1 Þ ¼ ps1;1 x2s1;3 " # F1 ðT1  xs1;2 Þ UAWm ðTj1  xs1;2 Þ þ þ ps1;2 : x2s1;3 dCp x2s1;3 ðA:31Þ

At the autonomous transition from ls1 to ls2 ; x ¼ xs1 ðts1 Þ ¼ xs2 ð0Þ ¼ ½xs1;1 xs1;2 xs1;3 T and p ¼ ps1 ðts1 Þ ¼ ½ps1;1 ps1;2 ps1;3 Tjt¼ts1 aps2 ð0Þ ¼ ½ps2;1 ps2;2 ps2;3 Tjt¼0 xs1;1 ¼ 2000:

p’s2;1 ¼

ðA:32Þ

3 qðxs1;1  2000Þ 7 6 qxs1;1 7 6 6 qðx  2000Þ 7 s1;1 7 6 ps2 ð0Þ ¼ ps1 ðts1 Þ ¼ 6 7ps1 7 6 qxs1;2 7 6 4 qðxs1;1  2000Þ 5 qxs1;3 2 3 ps1;1 ðts1 Þ  ps1 6 7 ¼ 4 ps1;2 ðts1 Þ 5 ps1;3 ðts1 Þ 2

p’s2;3 ¼

F2 ðC2  xs2;1 Þ ps2;1 x2s2;3 " # F2 ðT2  xs2;2 Þ UAWm ðTj2  xs2;2 Þ ps2;2 : þ þ x2s2;3 dCp x2s2;3 ðA:40Þ

At the controlled transition from ls2 to ls3 ; x ¼ xs2 ðts2 Þ ¼ xs3 ð0Þ ¼ ½xts2;1 xts2;2 xts2;3 T and p ¼ ps2 ðts2 Þ ¼ ps3 ð0Þ ¼ ½pts2;1 pts2;2 pts2;3 T Hs3 ð0Þ  Hs2 ðts1 Þ   F1 ðC1  xts2;1 Þ  F2 ðC2  xts2;1 Þ ¼ pts2;1 xts2;3  F1 ðT1  xts2;2 Þ  F2 ðT2  xts2;2 Þ þ pts2;2 xts2;3  UAWm ðTj1  Tj2 Þ þ dCp xts2;3 þ pts2;3 ðF1  F2 Þ ¼ 0:

ðA:41Þ

0ptpts3 ; x ¼ xs3 ¼ During phase ls3 ¼ l1 ; ½xs3;1 xs3;2 xs3;3 T and p ¼ ps3 ¼ ½ps3;1 ps3;2 ps3;3 T : x’ s3;1 ¼

F1 ðC1  xs3;1 Þ xs3;3

ðA:42Þ

x’ s3;2 ¼

F1 ðT1  xs3;2 Þ UAWm ðTj1  xs3;2 Þ þ xs3;3 dCp xs3;3

ðA:43Þ

gives the following equation: ps2;1 ð0Þ ¼ ps1;1 ðts1 Þ  ps1 :

ðA:33Þ

H2 ¼ H1 þ ðqðxs1;1  2000Þ=qtÞps1 gives the following equation:   F2 ðC2  xs1;1 Þ ðps2;1 ð0ÞÞ xs1;3   F1 ðC1  xs1;1 Þ  ðps1;1 þ ps1 Þ xs1;3  F2 ðT2  xs1;2 Þ  F1 ðT1  xs1;2 Þ þ ps1;2 xs1;3  UAWm ðTj2  Tj1 Þ þ dCp xs1;3 þ ps1;3 ðF2  F1 Þ ¼ 0:

ðA:34Þ

0ptpts2 ; x ¼ xs2 ¼ During phase ls2 ¼ l2 ; ½xs2;1 xs2;2 xs2;3 T and p ¼ ps2 ¼ ½ps2;1 ps2;2 ps2;3 T : F2 ðC2  xs2;1 Þ x’ s2;1 ¼ ðA:35Þ xs2;3 x’ s2;2 ¼

F2 ðT2  xs2;2 Þ UAWm ðTj2  xs2;2 Þ þ xs2;3 dCp xs2;3

x’ s2;3 ¼ F2  F

x’ s3;3 ¼ F1  F p’s3;1 ¼ p’s3;2 ¼

p’s3;3

F1 ps3;1 xs3;3  F1 UAWm þ ps3;2 xs3;3 dCp xs3;3

ðA:44Þ ðA:45Þ



ðA:46Þ

  F1 C1  xs3;1 ¼ ps3;1 x2s3;3 "    # F1 T1  xs3;2 UAWm Tj1  xs3;2 þ þ ps3;2 : x2s3;3 dCp x2s3;3 ðA:47Þ

in phase ls3 ¼ l1 ; x ¼ xs3 ðts3 Þ ¼ At t ¼ ts3 ½xs3;1 xs3;2 xs3;3 T and p ¼ ps3 ðts3 Þ ¼ ½ps3;1 ps3;2 ps3;3 T

ðA:36Þ

xs3;1 ¼ 2300:

ðA:48Þ

ðA:37Þ

xs3;2 ¼ 311:

ðA:49Þ

P. Manon et al. / Control Engineering Practice 10 (2002) 133–149

3 qðxs3;1  2300Þ qðxs3;2  311Þ 7 6 qxs3;1 qxs3;1 7" 6 # 6 qðx  2300Þ qðx  311Þ 7 p s3;1 s3;2 7 s3;1 6 ¼6 7 7 ps3;2 6 qxs3;2 qxs3;2 7 6 4 qðxs3;1  2300Þ qðxs3;2  311Þ 5 qxs3;3 qxs3;3 2 3 ps3;1 6 7 ¼ 4 ps3;2 5 2

ps3

0 gives the following equation: ps3;3 ¼ 0:

ðA:50Þ

References Asarin, A., Maler, O., & Pnueli, A. (1995). On the analysis of dynamical systems with piecewise constant derivatives. Theoretical Computer Science, 138, 35–65. Athans, M., & Falb, P. L. (1966). Optimal control: An introduction to the theory and its applications. New York: McGraw-Hill. Bonnard, B., & Pelletier, M. (1995). Time minimal synthesis for planar systems in the neighborhood of a terminal manifold of codimension one. Journal of Mathematical Systems, Estimation and Control, 5(3), 379–381. Branicky, M., Borkar, V., & Mitter, S. (1994). A unified framework for hybrid control. Proceedings of the 11th international conference on analysis and optimization of systems (pp. 352–358). SophiaAntipolis, France: Springer.

149

Bryson, A. E., & Ho, Y. C. (1969). Applied optimal control. Boston: Gin and Co.. ! Ciarlet, P. G. (1990). Introduction a" l’analyse numerique matricielle et a" l’optimisation, MASSON, Paris, 1990. ! es, ! Lavoisier (Chapter 19). Corriou, J. P. (1996). Commande des proced Culioli, J. C. (1994) Introduction a" l’optimisation, Ellipses. Johansson, K. H., Egerstedt, M., Lygeros, J., & Sastry, S. (1999). On the regularization of Zeno hybrid automaton. System and Control Letters, Special issue on hybrid systems, 38, 141–150. Manon, P. (2001). Sur l’optimisation des s!equences de fonctionement des syst"emes dynamiques hybrides, Ph.D. Thesis, University Claude Bernard, Lyon, France. Manon, P., & Valentin-Roubinet, C. (1999a). Synthesis of a trajectory to nominal mode for a class of hybrid dynamical systems. European Journal of Automation, 33(8–9), 995–1014. Manon, P., & Valentin-Roubinet, C. (1999b). Controller synthesis for hybrid systems with linear vector fields. Proceeding of the IEEE international symposium on intelligent control/intelligent systems and semiotics (ISIC 99) (pp. 17–22). Cambridge, US. Manon, P., Valentin-Roubinet, C., & Gilles, G. (2000). Optimal control of hybrid systems with the maximum principle: application to a non-linear process. IEEE conference on decision and control, Sydney, Australia. Pontryagin, L. S., Boltyanskii, V. G., Gamkrelidze, R. V., & Mishchenko, E. F. (1964). The mathematical theory of optimal processes. New York: Pergamon. Riedinger, P. (2000). Contribution a" la commande optimale des syst"emes dynamiques hybrides, Ph.D. Thesis, Centre de Recherche en Automatique de Nancy (CRAN), INPL, France. Riedinger, P., Zanne, C., & Kratz, F. (1999). Time Optimal Control of Hybrid Systems, proceedings of the ACC 99, San Diego, US. Tittus, M., & Egardt, B. (1998). Control design for integrator hybrid systems. IEEE Transactions on Automatic Control, 43(4), 491–500.