Exploitation of the control switching structure in multi-stage optimal control problems by adaptive shooting methods

Accepted Manuscript Title: Exploitation of the control switching structure in multi-stage optimal control problems by adaptive shooting methods Author...

Download PDF

601KB Sizes 0 Downloads 28 Views

Report

PDF Reader
Full Text

Accepted Manuscript Title: Exploitation of the control switching structure in multi-stage optimal control problems by adaptive shooting methods Author: Fady Assassa Wolfgang Marquardt PII: DOI: Reference:

S0098-1354(14)00327-5 http://dx.doi.org/doi:10.1016/j.compchemeng.2014.11.009 CACE 5073

To appear in:

Computers and Chemical Engineering

Received date: Revised date: Accepted date:

27-3-2014 23-11-2014 28-11-2014

Please cite this article as: Fady Assassa, Wolfgang Marquardt, Exploitation of the control switching structure in multi-stage optimal control problems by adaptive shooting methods, Computers and Chemical Engineering (2014), http://dx.doi.org/10.1016/j.compchemeng.2014.11.009 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Aachener Verfahrenstechnik, RWTH Aachen University, Turmstr. 46,52064 Aachen, Germany b German Research School for Simulation Sciences GmbH, 52425 J¨ ulich, Germany

an

us

a

cr

Fady Assassaa,b , Wolfgang Marquardta,∗

ip t

Exploitation of the control switching structure in multi-stage optimal control problems by adaptive shooting methods

Abstract

M

The adaptive switching structure approach is generalized from single-stage problems and single-shooting to multi-stage problems and multiple-shooting.

d

This generalization is based on previous work on the exploitation and detection of the control switching structure and wavelet-based control grid

te

refinement for single-stage problems. Here, single-shooting is employed to

Ac ce p

transcribe the multi-stage optimal control problem (OCP) into a nonlinear programming problem. The proposed multi-stage formulation is also capable to represent the transcription of single-stage OCP stemming from multipleshooting. Thus, the previously reported adaptive multiple-shooting approach is extended by an adaptation of the switching structure. Finally, a new stopping criterion is introduced that measures the intermediate constraint violation at the optimal solution. The proposed adaptive switching structure detection is illustrated for a ∗

Corresponding author. Tel.: +49 241 8094668 Fax: +49 241 8092326 Email address: [email protected] (Wolfgang Marquardt)

Preprint submitted to Computers and Chemical Engineering

December 3, 2014

Page 1 of 72

multi-stage and a multiple-shooting problem using the Williams-Otto semi-

ip t

batch reactor. A solution of user-specified accuracy in the objective and the

Keywords:

cr

path-constraints can be obtained using only few decision variables.

Dynamic optimization, control grid refinement, multi-stage optimal control,

1. Introduction

an

1

us

structure detection, multiple-shooting, single-shooting

The optimal operation of chemical processes often requires the solution

3

of multi-stage control problems. The scientific and industrial interest in

4

multi-stage optimal control problems (OCP) has been increasing steadily in

5

process systems engineering (Barton and Pantelides, 1994; Vassiliadis et al.,

6

1994a,b; Avraam et al., 1998; Bonvin, 1998; Biegler and Grossmann, 2004).

7

Generally, multi-stage OCP can be divided into two classes. The first class

8

assumes that the number of stages the system goes through is unknown,

9

while the second class requires a priori knowledge of the number and type

11

12

13

14

15

16

d

te

Ac ce p

10

M

2

of stages (Avraam et al., 1998). The first class often consists of multiple stages, where each stages is described by a diﬀerent DAE system. The stages are connected by state transitions, which are triggered by some logical constraint. These systems are also referred to as non-smooth systems1 . Typically, the non-smooth systems of interest are continuous in the state but with discontinuous right-hand side. These systems are modeled using if-else-clauses in the modeling environments. Technical examples for non1

For a detailed classification of non-smooth systems the reader is referred to the work

of (Leine and Nijmeijer, 2004).

2

Page 2 of 72

smooth systems include the model of a reactive three-phase batch distilla-

18

tion column by Br¨ uggemann et al. (2004) or just the model of a bouncing

19

ball (Hannemann-Tam´as et al., 2014).

ip t

17

Typical examples of the second class of multi-stage problems include the

21

operation of batch and semi-batch reactors, which rely on recipe-driven op-

22

eration (Bonvin, 1998). Recipes are subdivided hierarchically into process

23

stages, process operations and process actions (ISA, 2010). These opera-

24

tional modes usually result in a planned sequence of chemical or physical

25

changes of the state of material being processed. Thus, a batch process nor-

26

mally comprises various stages, e.g., a batch reactor is first cleaned, then

27

charged, the reaction is enabled, and the reactor is finally drained. De-

28

pending on the process, the reaction can be further subdivided into ad-

29

ditional stages, e.g., adding catalyst and reactants separately, withdraw-

30

ing material, heating and cooling (Bonvin, 1998). Thus, industrial semi-

31

batch processes always result in multi-stage OCPs with a known number

32

of stages. An example of particular industrial interest is the optimal con-

38

Ac ce p

te

d

M

an

us

cr

20

39

tion, in a dynamic cycle (Nilchan and Pantelides, 1998)). For these ap-

40

plications, each stage is described by separate state variables and mod-

41

els, which are linked by incorporating transitions between those dynamic

33

34

35

36

37

trol of (fed-batch) bioreactors (Balsa-Canto et al., 2000; Ramkrishna, 2003), which are used to manufacture specialty chemicals, pharmaceuticals, and bio-chemicals. Even continuous processes can experience transient operational stages (Barton and Pantelides, 1994) such as multiple grade changes. Further, chemical processes may include diﬀerent modes of operation requiring diﬀerent process models (e.g., adsorption, regeneration, pressuriza-

3

Page 3 of 72

stages (Biegler and Grossmann, 2004). A multi-stage OCP may also arise

43

in a control context, e.g., if an infinite horizon problem is formulated

44

(W¨ urth and Marquardt, 2014).

ip t

42

This paper deals with the formulation and optimization of multi-stage

46

OCPs arising in process engineering with a known number of stages. The

47

case of state transition networks is not discussed in this work, however, the

48

algorithm presented could be also generalized for this particular case. The

49

formulation of a multi-stage OCP with a given number of stages requires

50

setting up a sequence of operational modes (e.g. process stages, process

51

operations and process actions) with a varying set of controls, objectives,

52

constraints, and even model equations with complicated state transitions be-

53

tween the stages. Here, the stages of a multi-stage OCP are defined as either

54

physical stages, which are either enforced by state transition networks of the

55

underlying nonlinear process, or by algorithmic stages, which are due to the

56

formulation of the OCP (Schlegel and Marquardt, 2006a). To the author’s

57

knowledge, multi-stage OCP in the literature consists of either physical or

63

Ac ce p

te

d

M

an

us

cr

45

64

relying on a pragmatically chosen grid. However, it is important to note

65

that in such an approach the computed solution may not properly resolve

66

the qualitative structure of the optimal control profile (switching structure),

58

59

60

61

62

numerical stages. In this work, a generalized formulation that simultaneously comprises both physical and numerical stages is presented which results in a nested multi-stage problem. Typically, multi-stage OCPs are solved after approximating the contin-

uous problem by a finite-dimensional optimization problem by some kind of transcription. The control profiles are typically parameterized uniformly

4

Page 4 of 72

which is determined by a series of (continuous) arcs delimited by discon-

68

tinuous jumps (Bryson and Ho, 1975). Localizing the switching points be-

69

tween continuous arcs exactly not only improves the solution accuracy but,

70

more importantly, increases process insight. For example, this insight al-

71

lows applying NCO (necessary conditions of optimality) tracking control

72

(Srinivasan et al., 2003a) only using the solution structure determined by

73

oﬀ-line optimization (Kadam et al., 2007). A straightforward approach al-

74

lows to optimize the location of the grid points for a fixed number of control

75

intervals. Cuthrell and Biegler (1987) and von Stryk (1995) have applied this

76

idea in the context of full discretization methods. Vassiliadis et al. (1994a)

77

have used a similar concept, but apply it in the context of single shooting.

78

The major drawback of such relocation methods is that the resulting nonlin-

79

ear programming problem (NLP) tends to be strongly nonlinear and tough to

80

solve. In particular, even a linear optimal control problem results in an NLP

81

after parameterization if relocation methods are used. This drawback can be

82

compensated by introducing a bi-level approach, where the optimization is

88

Ac ce p

te

d

M

an

us

cr

ip t

67

89

multi-stage problems.

83

84

85

86

87

solved on a fixed grid in an inner loop, while the lengths of the discretization intervals are determined in an outer loop (Tanartkit and Biegler, 1997). To the authors’ knowledge, multi-stage OCPs cannot be tackled directly

by commercial software packages (e.g. gPROMS2 , Aspen Custom Modeler3 ). In this paper, we generalize the adaptive switching structure ap-

proach introduced by Schlegel and Marquardt (2006b) for single-stage to

2 3

This generalization is based on previous work

http://www.psenterprise.org/ http://www.aspentech.org/

5

Page 5 of 72

on the exploitation and detection of the control switching structure

91

(Schlegel and Marquardt, 2006a) and wavelet-based control grid refinement

92

(Schlegel et al., 2005) for single-stage problems. Here, we employ direct

93

single-shooting (Sargent and Sullivan, 1978) to transcribe the continuous dy-

94

namic multi-stage OCP into an NLP. Further, the proposed multi-stage for-

95

mulation is capable to represent the transcription of a single-stage OCP into

96

an NLP by multiple-shooting4 (Morrison et al., 1962; Bock and Plitt, 1984).

97

Thus, we also extend the adaptive multiple-shooting approach introduced by

98

Assassa and Marquardt (2014) by the adaptive switching structure approach

99

in this paper. Further, we introduce a new measure to evaluate the interme-

100

diate constraint violation (ICV), which can be used as a stopping criterion

101

for both refinement approaches.

M

an

us

cr

ip t

90

The adaptive switching structure approach of Schlegel and Marquardt

103

(2006b) generalized to multi-stage and multiple-shooting problems com-

104

prises two steps:

105

ated using the wavelet-based control grid refinement (Schlegel et al., 2005;

te

d

102

111

Ac ce p

First, a problem-dependent discretization is gener-

112

particular, every physical stage is subdivided by some algorithmic stages.

106

107

108

109

110

Assassa and Marquardt, 2014) to avoid an unnecessarily fine control grid. The iterative refinement procedure is terminated if the detected control switching structure for all stages remains unchanged in two consecutive iterations. In the second step, the multi-stage problem is reformulated to incorporate the detected switching structure (Schlegel and Marquardt, 2006a) by introducing a second layer of stages, the so-called algorithmic stages. In

4

Multiple-shooting explicitly addresses OCP problems embedding unstable initial value

problems, which often cannot be solved successfully by single-shooting.

6

Page 6 of 72

Each algorithmic stage has the same type of controls, objective, constraints

114

and model as its associated physical stage. These algorithmic stages allow to

115

capture the control switching structure behavior with only few optimization

116

parameters combining control grid adaptation and switching structure de-

117

tection. The control grid of the reformulated problem is again refined using

118

wavelet-based grid refinement until a satisfactory solution quality is obtained.

119

The stopping criterion either relies on a measure of the ICV and/or a change

120

in the objective function. A change in the switching structure leads to a

121

reformulation of the underlying problem with subsequent control grid refine-

122

ment.

an

us

cr

ip t

113

We illustrate adaptive control grid refinement and adaptive switching

124

structure detection for a multi-stage problem and a multiple-shooting strat-

125

egy using the Williams-Otto semi-batch reactor (Williams and Otto, 1960).

126

The proposed algorithm successfully detects the control switching structure

127

in both cases. A solution of user-specified accuracy in the objective and the

128

path-constraints can be obtained using only few decision variables.

134

Ac ce p

te

d

M

123

135

strategy. Section 5 introduces an exemplary multi-stage and single-stage

136

problem of the Williams-Otto semi-batch reactor. The multi-stage problem

137

is solved using single-shooting, while the single-stage problem is discretized

129

130

131

132

133

The paper is structured as follows. Section 2 first introduces the multi-

stage problem formulation and the necessary conditions of optimality. The theoretical background of the suggested approach and the applied solution method are presented in Section 3. It concludes with a summary of the unified automatic structure detection approach. In Section 4, we give an overview on the solvers used and the implementation of the optimization

7

Page 7 of 72

by the multiple-shooting strategy. Both examples are analyzed and compared

139

to an equidistant discretization of the controls with respect to performance

140

and solution accuracy. Finally, in Section 6, we conclude the paper with a

141

summary and future perspectives.

142

2. Problem formulation

us

cr

ip t

138

First, the multi-stage OCP is formulated and generalized by introducing

144

a second level of stages, so-called algorithmic stages. Second, the theoretical

145

framework to reformulate the multi-stage OCP to incorporate the solution

146

structure using algorithmic stages is introduced.

147

2.1. Multi-stage optimal control problem

M

an

143

The general multi-stage formulation of a dynamic optimization (or opti-

149

mal control) problem allowing for a changing set of model equations, controls,

150

constraints and objectives from one stage to another reads as follows:

Ac ce p

te

d

148

min

Φ=

uk (t),∆tk ,xk,0 ,∀k∈K

K ∑

Φk (xk (tk ), zk (tk ))

(1a)

k=1

s.t. Mk x˙ k (t) = fk (xk (t), zk (t), uk (t)),

∀k ∈ K,

(1b)

0 = gk (xk (t), zk (t), uk (t)),

∀k ∈ K,

(1c)

0 = xk (tk−1 ) − xk,0 ,

∀k ∈ K,

(1d)

∀k ∈ K,

(1e)

∀k ∈ K\K,

(1f)

0 ≤ hk (xk (t), zk (t)),

∀k ∈ K,

(1g)

0 ≤ ek (xk (tk ), zk (tk )),

∀k ∈ K,

(1h)

Up xLo k,0 ≤ xk,0 ≤ xk,0 ,

0 = xk (tk ) − xk+1 (tk ),

8

Page 8 of 72

∀k ∈ K,

(1i)

tk = tk−1 + ∆tk ,

∀k ∈ K,

(1j)

∀k ∈ K,

(1k)

ip t

Up uLo k ≤ uk (t) ≤ uk ,

Up ∆tLo k ≤ ∆tk ≤ ∆tk ,

≤

K ∑

∆tk ≤ ∆tU p ,

cr

∆t

Lo

k=1

∀k ∈ K.

us

t ∈ [tk−1 , tk ],

(1l) (1m)

K is the number of (physical) stages and the set K = {1, . . . , K} com-

152

prises the indices of all stages. Thus, the considered time horizon [0, tf ] is

153

partitioned such that 0 = t0 < t1 < · · · < tk < · · · < tK = tf . The

154

DAE model (Eqs. (1b)-(1c)) describing the system behavior on the k-th

155

stage is represented by the nonlinear functions fk : Rnx × Rnz × Rnu → Rnx ,

156

gk : Rnx ×Rnz ×Rnu → Rnz and the non-singular mass matrix Mk . We assume

157

Mk to be constant and the DAE systems to have a diﬀerential index of at

158

most 1. Note that high-index DAE systems can be reduced to index-one DAE

159

systems by one of the existing symbolical algorithms (e.g. Pantelides, 1988;

k

k

k

k

165

Ac ce p

te

d

k

M

an

151

166

(1f). xk,0 refers to the vector of consistent initial conditions at time t = tk−1 .

167

If the initial values xk,0 , ∀k ∈ K\1 are decision variables of the OCP, then

168

we can interpret the (physical) stage k ∈ K as a shooting interval required to

160

161

162

163

164

Bachmann et al., 1990; Mattsson and S¨oderlind, 1993; Unger et al., 1995). Here, xk (t) ∈ Rnx and zk (t) ∈ Rnz denote the vectors of diﬀerential and k

algebraic states of stage k ∈ K , respectively. For the purpose of this paper, it is assumed that the number of diﬀerential variables and equations is identical for all stages k ∈ K . Thus, continuous diﬀerential state variables are enforced across stage boundaries by the simple stage transition condition

9

Page 9 of 72

formulate a multiple-shooting problem. However, if these initial values are

170

not decision variables, then the initial condition (1d) is only valid for stage

171

k = 1 resulting in a transcription of single-shooting type.

ip t

169

172

In addition to the initial values xk,0 , the time-dependent control variables

173

uk (t) ∈ Rnu and the stage lengths ∆tk are the decision variables of the

174

continuous multi-stage optimization problem. Without loss of generality,

175

the objective functions Φk : Rnx × Rnz → R are of Mayer-type, since an

176

integral cost can be reformulated as a terminal cost by including additional

177

diﬀerential variables and equations to the DAE system (Hazewinkel, 2001, cf.

178

Bolza problem). Note that it is not required to formulate a non-zero objective

179

function for each stage, however, at least one non-vanishing objective function

180

has to be given. Furthermore, state path-constraints (1g) are defined on the

181

stage horizon ∆tk whereas the terminal constraints (1h) are only formulated

182

at final time tk of stage k ∈ K . Simple bounds on the decision variables are

183

given in Eqs. (1e), (1i), and (1k), respectively, and a linear constraint on the

184

time horizon of the entire problem is given in Eq. (1l). Note that combined

186

187

188

189

190

M

d

te

Ac ce p

185

an

k

us

cr

k

state and control constraints can be converted to the considered formulation by including additional variables and equations to the DAE system (e.g., Hannemann-Tam´as (2012)). For brevity and without loss of generality, we do not consider time-invariant system parameters, which could easily be added and partly included in the set of decision variables. 2.2. Necessary conditions of optimality and switching structure

191

The necessary conditions of optimality (NCO) for the multi-stage OCP

192

can be derived from Pontryagin’s Minimum Principle (PMP). However,

193

for the sake of brevity, we omit deriving the NCO and refer the reader 10

Page 10 of 72

to Bryson and Ho (1975). Here, we only give a qualitative interpretation

195

of the NCO and their eﬀects on the solution structure as described by

196

Srinivasan et al. (2003b) and Schlegel and Marquardt (2006a).

ip t

194

cr

The course of an optimal control profile can be divided into four types of

197

diﬀerent arcs:

199

uUi,kp : the control is at its upper bound;

200

uLo i,k : the control is at its lower bound;

201

uPi,kath : the control is determined by an active path-constraint;

202

uSens : the control is sensitivity-seeking and thus minimizes the objective i,k

M

an

us

198

function.

203

The optimal solution structure of a multi-stage OCP always consists of a

205

distinct sequence of these arcs, which forms the so-called switching struc-

206

ture.

207

tact or touch-and-go points, where a constraint is active only in isolated

209

210

211

212

213

214

te

Additionally, the control might be determined by so-called con-

Ac ce p

208

d

204

points (Jacobson et al., 1971). Arcs and touch-and-go points exist depending on the order of the corresponding state path-constraint5 (Bryson and Ho, 1975). A contact point may be created for state path-constraints of arbitrary order by choosing the corresponding bounds in a specific way (Maurer and Heidemann, 1975). If such a bound is further tightened, a contact point transforms into an arc with active state path-constraint. Touchand-go points may transform into arcs for state path-constraints of even 5

A pure state path-constraint has order q, if an explicit expression for a control can be

obtained at the q-th time diﬀerentiation of the corresponding state path-constraint.

11

Page 11 of 72

Table 1: Touch-and-go points and arcs for diﬀerent order q of the state path-constraint

q-th time derivative touch-and-go points arcs

no

2,4,6,8,...

yes

3,5,7,9,...

yes

cr

1

yes yes

us

no

yes no

an

0

ip t

(Jacobson et al., 1971; B¨ uskens, 1998).

order (q ≥ 2) but remain touch-and-go points for odd-ordered state path-

216

constraints (q ≥ 3) regardless of the chosen constraint (Speyer et al., 1971).

217

An overview of the occurrence of touch-and-go points and arcs for OCP is

218

given in Table 1 (B¨ uskens, 1998). The analysis of the switching structure is

219

limited to arcs only corresponding to an order of the state path-constraints of

220

at most one, since touch-and-go points rarely occur in practical application

221

(Schlegel and Marquardt, 2006a).

223

224

225

226

227

228

d

te

Ac ce p

222

M

215

Schlegel and Marquardt (2006a) automatically estimate the optimal

switching structure without requiring information on the Hamiltonian function or the adjoint system. Their analysis carries over to multi-stage OCPs to allow for an automated reformulation to incorporate the optimal switching structure on each of the physical stages as explained in the following subsection.

2.3. Algorithmic stages

229

In this subsection, we reformulate the multi-stage OCP (1) using a second

230

layer of stages, which is embedded into all physical stages k ∈ K . These 12

Page 12 of 72

underlying algorithmic stages transcribe the OCP to cope with the optimal

232

switching structure. The OCP obtained by embedding these algorithmic

233

stages introduces additional linear constraints and is as follows:

cr

Φ=

Φk (xLk k (tk ), zLk k (tk ))

k=1

us

min

ulk (t),∆τkl ,xk,0

K ∑

ip t

231

(2a)

∀k ∈ K,

∀l ∈ Lk ,

(2b)

0 = gk (xlk (τ ), zlk (τ ), ulk (τ )),

∀k ∈ K,

∀l ∈ Lk ,

(2c)

an

s.t. Mk x˙ lk = fk (xlk (τ ), zlk (τ ), ulk (τ )),

0 = x1k (tk−1 ) − xk,0 ,

(2d)

∀k ∈ K,

(2e)

∀k ∈ K\K,

(2f)

M

Up xLo k,0 ≤ xk,0 ≤ xk,0 ,

∀k ∈ K,

∀k ∈ K,

∀l ∈ Lk \Lk ,

(2g)

d

0 = xLk k (tk ) − x0k+1 (tk ), l 0 = xlk (τkl ) − xl+1 k (τk ),

∀k ∈ K,

∀l ∈ Lk ,

(2h)

0 ≤ ek (xLk k (tk ), zLk k (tk )),

∀k ∈ K,

te

0 ≤ hk (xlk (τ ), zlk (τ )),

Ac ce p

Up l uLo k ≤ uk (τ ) ≤ uk ,

∆τkLo ≤ ∆τkl ≤ ∆τkU p = ∆tUk p , τkl = τkl−1 + ∆τkl ,

∆tk =

Lk ∑

∆τkl ,

(2i)

∀k ∈ K,

∀l ∈ Lk ,

(2j)

∀k ∈ K,

∀l ∈ Lk ,

(2k)

∀k ∈ K,

∀l ∈ Lk ,

(2l)

∀k ∈ K,

(2m)

∀k ∈ K,

(2n)

∀k ∈ K,

(2o)

l=1

tk = tk−1 + ∆tk ,

Up ∆tLo k ≤ ∆tk ≤ ∆tk ,

∆t

Lo

≤

K ∑

∆tk ≤ ∆tU p ,

(2p)

k=1

∀k ∈ K,

τ ∈ [τkl−1 , τkl ],

∀l ∈ Lk

(2q)

13

Page 13 of 72

We have introduced Lk ≥ 1 algorithmic stages for each physical stage k ∈ K.

235

The continuity of the diﬀerential state variables across algorithmic stages

236

is ensured by Eq. (2g) in addition to the continuity condition across stage

237

boundaries, Eq. (2f). The reformulated problem introduces the length of the

238

algorithmic stage ∆τkl as decision variables to the OCP. We have to choose

239

the lower bound of the algorithmic stage (cf. Eq. (2k)) such that τkLo > 0

240

to avoid numerical instability (Sager, 2009). The stage durations ∆tk are

241

determined by the equality constraint (2m).

242

3. Numerical solution

M

an

us

cr

ip t

234

The solution strategies for OCP can be classified into indirect and di-

244

rect methods. Binder et al. (2001a) and Srinivasan et al. (2003b) have re-

245

viewed the specific advantages and disadvantages of the two approaches in

246

detail. Here, we only consider direct methods since indirect methods are

247

diﬃcult to apply to OCPs encountered in chemical engineering applications

248

(Binder et al., 2001a; Biegler, 2010).

250

251

252

253

254

te

Ac ce p

249

d

243

Direct methods approximate the infinite-dimensional OCP by a finite-

dimensional (algebraic) optimization problem, which is solved by means of a suitable NLP algorithm. Three diﬀerent types of direct approaches have been established: single-shooting, also referred to as control-vector parameterization (Sargent and Sullivan, 1978), multiple-shooting (Morrison et al., 1962; Bock and Plitt, 1984), and simultaneous state- and control-vector pa-

255

rameterization, often called full discretization approach (Tsang et al., 1975;

256

Biegler, 1984; Kraft, 1985; Biegler, 2010). These approaches diﬀer in the way

257

the variables of the OCP are discretized. 14

Page 14 of 72

In this work, direct single- and multiple-shooting are employed to solve

259

the OCP. Either of the two solution strategies can be applied to the general

260

problem formulation introduced in Eq. (2) by either ignoring or consider-

261

ing the initial conditions xk,0 , ∀k ∈ K\1, as decision variables. In both

262

approaches, the control variables ulk (t) have to be parameterized likewise.

263

3.1. Parameterization of control variables and state path-constraints

us

cr

ip t

258

In single- and multiple-shooting the continuous control variables ulk (t)

265

are explicitly parameterized. This can, for example, be done using B-splines

266

(de Boor, 1978) of order r, such that the parameterized controls u˜lk,i with

267

i = 1, . . . , nku are represented by

j∈∆ul

269

270

i = 1, . . . , nku ,

∀k ∈ K, ∀l ∈ Lk ,

(3)

with decision parameters clk,i,j parameterizing the nku controls ulk (t) by a finite [ ] set ∆ulk,i ⊂ τkl−1 , τkl of time-dependent basis functions φl,r k,j (t). B-splines of

Ac ce p

268

clk,i,j φl,r k,j (τ ),

te

k,i

M

∑

d

u˜lk,i (clk,i , τ ) =

an

264

order r are given by

φl,1 k,j (τ ) :=

φl,r k,j (τ ) :=

  1 if ξj ≤ τ ≤ ξj+1  0 else

j ∈ ∆ulk,i ,

(4)

τ − ξj ξj+r − τ l,r−1 φl,r−1 φ (τ ), r > 1, j ∈ ∆ulk,i k,j (τ ) + ξj+r−1 − ξj ξj+r − ξj+1 k,j+1 (5)

271

and are centered on a chosen set of knots ξj . In this work, the approximation

272

order r is chosen to be 1 (piecewise constant) or 2 (piecewise linear) to prevent 15

Page 15 of 72

violation of path-constraints (Schlegel et al., 2005). In the following, the

274

control parameters clk,i,j are collected in the vector

275

clk = [clk,1,1 , . . . ,clk,1,j ,. . . ,clk,nk ,j ,. . . clk,nk ,j ,. . . ]T . u

The state path-constraints (2h) are discretized on a finite set ∆hlk ⊂ ] τkl−1 , τkl as

us

277

[

u

cr

276

ip t

273

278

an

0 ≤ h(xlk (τj ; clk ), z(τj ; clk )), j ∈ ∆hlk ,

(6)

which is chosen twice as fine as the control parameterization grid ∆ulk,i ⊂ ∆lhl

k

(Schwartz, 1996). Though this strategy does not guarantee intermediate

280

path-constraint violations to disappear, it still often reduces them in magni-

281

tude as it will be demonstrated in Section 5.

282

3.2. NLP resulting from single-shooting

d

M

279

Transcribing the continuous OCP into an NLP using single-shooting as-

284

sumes a continuous state trajectory for all (physical) stages k ∈ K . Thus,

286

287

288

289

290

Ac ce p

285

te

283

the state continuity condition (2f) and (2g) are implicitly fulfilled by the numerical integration routine for all stage transitions. Consequently, we can discard the initial conditions xk,0 for all stages k ∈ K\1, if the initial conditions x1,0 are given. Hence, the decision variables consist of the control parameters clk and the algorithmic stage durations ∆τkl which are concatenated to plk = [(clk )T , (∆τkl )T ]T . Then, the resulting NLP can be written as

291

min Φ = plk

K ∑

Φk (xLk k (tk ; plk ), zLk k (tk ; plk ))

(7a)

k=1

16

Page 16 of 72

s.t. 0 ≤ hk (xlk (τj ; plk ), zlk (τj ; plk )),

∀k ∈ K,

∀l ∈ Lk ,

∀k ∈ K,

Up l uLo k ≤ ck ≤ uk ,

∀k ∈ K, ∀k ∈ K, ∀k ∈ K,

τkl = τkl−1 + ∆τkl , ∆τkl ,

l=1

Up ∆tLo k ≤ ∆tk ≤ ∆tk ,

k=1

∀l ∈ Lk ,

(7e)

∆tk ≤ ∆tU p .

∀l ∈ Lk ,

(7f) (7g)

∀k ∈ K,

(7h)

∀k ∈ K,

(7i) (7j)

d

K ∑

M

tk = tk−1 + ∆tk ,

∆tLo ≤

(7d)

∀k ∈ K,

an

∆tk =

Lk ∑

∀l ∈ Lk ,

us

∆τkLo ≤ ∆τkl ≤ ∆τkU p = ∆tUk p ,

(7c)

cr

0 ≤ ek (xLk k (tk ; plk ), zLk k (tk ; plk )),

(7b)

ip t

j ∈ ∆lh,k ,

Note that the bounds on the control (Eq. (7d)) are directly enforced for

293

the control parameters clk . This is only valid for piecewise constant and/or

294

piecewise linear control approximations, which are exclusively used in this

299

Ac ce p

te

292

300

intervals. Their number is typically chosen such that numerical stability is

301

ensured (Bock, 1987; Hartwich, 2013) and that the computational load of

302

the NLP is decreased (Mattheij and Staarink, 1984). However, algorithmic

295

296

297

298

paper.

3.3. NLP resulting from multiple-shooting After control vector parameterization, we can interpret the general multi-

stage formulation Eq. (1) also as a multiple-shooting discretization of a single-stage OCP. The physical stages K become equivalent to the shooting

17

Page 17 of 72

stages are still used to cope with the characteristic switching structure of the

304

OCP.

ip t

303

Multiple-shooting requires free initial conditions xk,0 for stages k ∈ K\1,

306

assuming known initial conditions x1,0 . Therefore, we have to include the

307

boundary condition (2f) into the NLP to enforce continuity across physical

308

stage or shooting interval boundaries tk with k ∈ K\1. The continuity condi-

309

tion across algorithmic stage boundaries (2g) is again enforced by numerical

310

integration.

an

us

cr

305

The decision variables for the multiple-shooting formulation consist of the

312

control parameters clk , the algorithmic stage durations ∆τkl and the initial

313

conditions xk,0 for k ∈ K\1 which are concatenated to

314

plk = [(clk )T , (∆τkl )T , (xk,0 )T ]T . Then, the resulting NLP can be written as

plk ,∀k∈K

K ∑ k=1

Φk (xLk k (tk ; plk ), zLk k (tk ; plk ))

te

min Φ =

d

M

311

Ac ce p

s.t. 0 ≤ hk (xlk (τj ; plk ), zlk (τj ; plk )),

(8a) ∀k ∈ K,

∀l ∈ Lk , j ∈ ∆lh,k ,

0 = x1k (tk−1 ) − xk,0 ,

Lo xk,0 ≤ xk,0 ≤ xUk,0p

0 = xLk k (tk ) − x0k+1 (tk ), 0 ≤ ek (xLk k (tk ; plk ), zLk k (tk ; plk )),

Up l uLo k ≤ ck ≤ uk ,

∆τkLo ≤ ∆τkl ≤ ∆τkU p = ∆tUk p , τkl = τkl−1 + ∆τkl ,

(8b)

∀k ∈ K\1,

(8c)

∀k ∈ K\1,

(8d)

∀k ∈ K\K,

(8e)

∀k ∈ K,

(8f)

∀k ∈ K,

∀l ∈ Lk ,

(8g)

∀k ∈ K,

∀l ∈ Lk ,

(8h)

∀k ∈ K,

∀l ∈ Lk ,

(8i)

18

Page 18 of 72

∆tk =

Lk ∑

∀k ∈ K,

∆τkl ,

(8j)

∀k ∈ K,

tk = tk−1 + ∆tk ,

∆t

Lo

≤

K ∑

∀k ∈ K,

cr

Up ∆tLo k ≤ ∆tk ≤ ∆tk ,

∆tk ≤ ∆tU p .

(8k)

(8l)

(8m)

us

k=1 315

ip t

l=1

3.4. Wavelet-based adaptive control vector parameterization

A suﬃciently accurate control profile is essential to automatically detect

317

the switching structure of the optimal control profile and to obtain reason-

318

able approximations for sensitivity-seeking or path-constrained arcs. In our

319

previous publication (Assassa and Marquardt, 2014), we addressed how both

320

these problems can be tackled by automatic adaptive control grid refinement.

321

Thus, we limit our introduction to the method that we employ in this paper.

322

Here, we apply the wavelet6 -based adaptive refinement strategy intro-

323

duced by Binder et al. (2000) and extended by Schlegel et al. (2005) and

324

Assassa and Marquardt (2014) for single- and multiple-shooting, respec-

326

327

328

329

330

M

d

te

Ac ce p

325

an

316

tively. This strategy starts from an initially coarse parameterization and refines the control grid using a multi-scale analysis of the previously obtained optimal solution. Additional grid points are only inserted where required and redundant grid points are eliminated to meet user-defined tolerances. The B-spline basis used to discretize the control variable has the advan-

tage that it relates to a certain class of wavelet functions. These wavelets 6

A comprehensive description of wavelets and the theory of multi-scale methods is be-

yond the scope of this article. The reader is referred to the work of Binder et al. (2001b,c) and the excellent review of Dahmen (1997).

19

Page 19 of 72

are well-suited for a multi-scale analysis, which is required to employ the

332

wavelet-based refinement strategy. The multi-scale analysis allows describ-

333

ing a function by a hierarchical representation with a collection of (wavelet)

334

coeﬃcients. Thus, the local frequency information of a function is eﬃciently

335

captured by these coeﬃcients, which are obtained by applying a wavelet

336

transformation to the discretized controls u˜lk,i (clk,i , τ ) (cf. Eq. (3)). The

337

wavelet coeﬃcients dlk,i obtained from the wavelet transformation are a mea-

338

sure for the local frequency content of a function. This allows for the for-

339

mulation of the necessary criteria to automatically adapt the control grid

340

using the wavelet coeﬃcients dlk,i as well as the equivalence between the

341

norm of the discretized controls u˜lk,i (clk,i , τ ) and the wavelet coeﬃcients dlk,i

342

(Schlegel et al., 2005). Small wavelet coeﬃcients can be discarded, because

343

they will only cause small changes in the approximate representation of the

344

original function ulk,i (τ ). Secondly, we can approximate functions, which

345

show strong local variations by linear combinations of relatively few wavelet

346

basis functions corresponding to the significant wavelet coeﬃcients.

348

349

350

351

cr

us

an

M

d

te

Ac ce p

347

ip t

331

The reader is referred to Schlegel et al. (2005) for a more detailed pre-

sentation of the algorithmic steps of the wavelet-based adaptive refinement procedure.

3.5. Automatic detection of the switching structure Schlegel and Marquardt (2006a, 2006b) describe an algorithm that auto-

352

matically detects the control switching structure, which typically appears in

353

the solution of OCPs (cf. Section 2.2). However, their method is limited to

354

OCPs consisting of a single physical stage. The key steps of the automatic

355

detection of the switching structure presented by Schlegel and Marquardt 20

Page 20 of 72

356

(2006a) are 1. the solution of the OCP by an appropriate discretization,

358

2. the detection of arcs in the optimal control profile using the NCO (cf.

361

cr

360

Section 2.2), and

3. the reparameterization as a multi-stage problem employing algorithmic stages using the optimal control trajectory from Step 2.

us

359

ip t

357

The advantage of reparameterization (in Step 3) is that a numerical solution

363

of higher accuracy and a better conditioning with fewer decision variables

364

can be obtained. Schlegel and Marquardt (2006a) achieve this by parame-

365

terizing the control variables with active control bounds (uU p , uLo ) by piece-

366

wise constant functions, which are excluded from the decision-variable vec-

367

tor. Control variables with inactive control bounds with sensitivity-seeking

368

or path-constrained arcs are approximated by piecewise linear, continuous

369

functions. Furthermore, Step 3 introduces a new algorithmic stage for every

370

detected arc, where the length of the algorithmic stages are additional deci-

372

M

d

te

Ac ce p

371

an

362

sion variables in the optimization. Thus, the original single-stage problem is transcribed into a multi-stage problem as follows:

min Φ = Φ1 (xL1 1 (t1 ), zL1 1 (t1 ))

(9a)

ul1 (τ ),∆τ1l

s.t. M1 x˙ l1 = f1 (xl1 (τ ), zl1 (τ ), ul1 (τ )),

∀l ∈ L1 ,

(9b)

0 = g1 (xl1 (τ ), zl1 (τ ), ul1 (τ )),

∀l ∈ L1 ,

(9c)

0 = x11 (t0 ) − x1,0 ,

(9d)

l 0 = xl1 (τ1l ) − xl+1 1 (τ1 ),

∀l ∈ L1 \L1 ,

(9e)

21

Page 21 of 72

0 ≤ h1 (xl1 (τ ), zl1 (τ )),

∀l ∈ L1 ,

0 ≤ e1 (xL1 1 (t1 ), zL1 1 (t1 )),

ip t

(9g)

Up l uLo 1 ≤ u1 (τ ) ≤ u1 ,

∆τ1l ,

l=1

an

t1 = t1−1 + ∆t1 ,

M

Up ∆tLo 1 ≤ ∆t1 ≤ ∆t1 ,

τ ∈ [τ1l−1 , τ1l ],

(9h)

∀l ∈ L1 ,

(9i)

∀l ∈ L1 ,

(9j)

us

τ1l = τ1l−1 + ∆τ1l , L1 ∑

∀l ∈ L1 ,

cr

∆τ1Lo ≤ ∆τ1l ≤ ∆τ1U p = ∆tU1 p ,

∆t1 =

(9f)

(9k) (9l) (9m)

∀l ∈ L1 .

(9n)

The multi-stage formulation (9) introduced by Schlegel and Marquardt

374

(2006a) corresponds to the previously introduced multi-stage formulation

375

with underlying algorithmic stages (2), if one physical stage (K = 1) and

376

L1 ≥ 1 algorithmic stages are introduced7 . This formulation allows for a

378

379

380

381

382

te

Ac ce p

377

d

373

precise detection of the switching structure of a single-stage problem if a suﬃciently accurate control profile is available (cf. Section 3.4). The automatic switching structure detection is illustrated for an exem-

plary solution trajectory in Figure 1 of a problem with a single physical stage (i.e., K = 1). First (Figure 1, top), the control intervals of the optimal solution are categorized according to the NCO (cf. Section 2.2). Neigh7

Note that it is easy to show that the automatic structure detection algorithm intro-

duced by Schlegel and Marquardt (2006a) can be formulated using K ≥ 1 physical stages and Lk = 1 algorithmic stages. However, then the interpretation used in this paper of physical and algorithmic stages becomes misleading.

22

Page 22 of 72

boring intervals of the same type are grouped together such that they form

384

an arc of the same type. Second (Figure 1, middle), an algorithmic stage

385

is introduced for every appearing arc resulting in Lk algorithmic stages.

386

The information about the type of arcs is stored in the Boolean vector

387

Ak,ui ∈ RLk and i = 1, . . . , nku . LK = 4 algorithmic stages are introduced with

388

A1,u1 = [Lo, Lo, U p, U p] and A1,u2 = [U p, P ath, P ath, Lo]. Finally (Figure 1,

389

bottom), the controls are reparameterized according to the detected switching

390

structure. Note that the decision variables only consist of the reparameter-

391

ized control with active path-constraints and (if present) sensitivity-seeking

392

arcs as well as the length of the newly introduced stages τkl , ∀l ∈ Lk . Thus,

393

we obtain a solution of a better-conditioned problem with fewer decision vari-

394

ables but higher accuracy in the subsequent optimization step. For technical

395

details on the algorithm, we refer to Schlegel and Marquardt (2006a,b).

d

M

an

us

cr

ip t

383

The extension of the automatic structure detection from single-stage to

397

multi-stage and multiple-shooting OCP is conceptually straightforward. The

398

main diﬃculty lies in providing an adequate implementation, which handles

404

Ac ce p

te

396

405

model equations with complicated state transitions between both stages.

406

This change in the problem formulation from one stage to the next typ-

407

ically leads to a change in the control switching structure which is han-

399

400

401

402

403

a multi-stage formulation with underlying algorithmic stages (2). We do not have to distinguish between multi-stage problems solved by single-shooting and single-stage problems solved by multiple-shooting using physical stages. First, for multi-stage problems solved by single-shooting, the transition

from one physical stage to another corresponds to a beginning of a new operational phase with a varying set of controls, objective, constraints, and even

23

Page 23 of 72

Up

Up

control u

0.8

2

Path

Lo

Lo 0.1

Lo 0.2

Lo 0.3

Lo

Lo

0.4

0.5 time

Lo 0.6

Lo 0.2

0.3

∆t

0.7

0.8

∆t

d

2

0.4

te

0.3

0.5 time

Lo 0.9

1

Up

∆t

3

4

Lo

0.2

1

Path

∆t

1

Lo 0.1

0.6

M

Path

0.4

0

0.5 time

Up

0.6

0

0.4

Up

0.8

0.9

Up

an

Lo 0.1

Lo

us

Path

0.2

0.2

0.8

Path

0.4

1

0.7 Up

0.6

0

Lo

Up

0.8

0

ip t

Path Path

0.2

1

Up

Path

0.4

0

Up

Path

Path

0.6

0

Up

control u

1

cr

1

Lo 0.6

0.7

0.8

0.9

1

Ac ce p

Figure 1: Automatic switching structure detection for a case with two controls u1 (t) and u2 (t).

408

409

410

411

dled by the underlying set of algorithmic stages. Thus, treating all physical stages k ∈ K independently by the structure detection algorithm of Schlegel and Marquardt (2006a) is reasonable. Second, physical stages interpreted as multiple-shooting intervals are de-

412

coupled due to the free initial values at the shooting nodes. Multiple-shooting

413

only enforces state continuity at the optimal solution using the stage bound-

414

ary condition (8c). The continuity condition in the state variables usually

415

leads to continuous controls at the stage boundaries. However, in rare cases 24

Page 24 of 72

the control variable may be discontinuous, if the stage boundary coincides

417

with a changing switching structure in the control variable. Distinguishing

418

between multiple-shooting intervals that have continuous and discontinuous

419

control variables at stage boundaries leads to an algorithm with increased

420

complexity, which is not universally applicable to multi-stage problems. An

421

algorithm that explicitly treats the continuity property of the control vari-

422

ables at stage boundaries will result in smaller NLP problems. Here, the

423

increase in the NLP size is accepted in return of a more general structure

424

detection approach. The load for sensitivity calculations introduced by the

425

additional decision variables is low due to the decoupling of multiple-shooting.

426

Further, our formulation allows to mix physical stages interpreted as

427

multiple-shooting intervals and physical stages implementing a multi-stage

428

formulation with altering objective function, constraints and controls. A

429

multiple shooting node is easily obtained by including all initial variables

430

xk,0 for some stage k ∈ K into the decision variable vector plk (cf. NLP

431

(8)). Algorithmic stages are used equivalently in both cases to cope with the

433

434

435

436

cr

us

an

M

d

te

Ac ce p

432

ip t

416

detected switching structure. 3.6. Unified automatic structure detection and control grid adaptation Here, we first present the unified automatic structure detection

and control grid adaptation for single-stage problems as introduced by Schlegel and Marquardt (2006b). Second, the algorithm for multi-stage and

437

multiple-shooting OCP is extended. The entire solution is summarized in

438

Algorithm 1 and 2 at the end of this subsection.

25

Page 25 of 72

439

3.6.1. Single-stage problems The algorithm of Schlegel and Marquardt (2006b) consists of mainly two

441

steps: The first step iteratively adapts the single-stage solution until the

442

detected switching structure Ak between two successive refinement steps re-

443

mains unchanged. Typically, the detected switching structure and the initial

444

guess of the switching times are well-suited to initialize a single-stage problem

445

with underlying algorithmic stages. A control grid discretization is suﬃcient

446

as long as the switching structure remains the same and the switching points

447

are suﬃciently close to the optimum (Sager, 2009).

an

us

cr

ip t

440

In the second step, the reformulated multi-stage problem is solved and

449

subsequently the initially coarsely discretized control grid of sensitivity-

450

seeking and path-constrained arcs is refined. Each sensitivity-seeking and

451

path-constrained arc is again analyzed for any structural changes, which may

452

result in the introduction of additional arcs. By construction, each algorith-

453

mic stage represents one arc (e.g., an initially detected sensitivity-seeking

454

arc), such that three diﬀerent cases are to be distinguished:

456

457

458

459

d

te

Ac ce p

455

M

448

1. The number of arcs remains the same, while having the same type. 2. The number of arcs remains the same, but some arc may change its type.

3. The number of arcs changes, which indicates a possible change of the switching structure.

460

The first case is typical for adaptively refined control grids. The second case

461

leads to the elimination of one of two neighboring algorithmic stages with

462

matching type. Finally, the third case requires a reformulation of the problem

463

to include the emerging arcs and their corresponding switching times. 26

Page 26 of 72

465

The unified automatic structure detection and control grid adaptation approach presented by Schlegel and Marquardt (2006b) stops, when

ip t

464

|Φκ − Φκ−1 | > ϵΦ , |Φκ |

cr

(10)

i.e. when the objective function Φ does not improve significantly during the

467

course of refinement steps κ. Here, ϵΦ is a user-specified stopping threshold.

468

3.6.2. Extension to multi-stage and multiple-shooting problems

us

466

The extension of the framework presented for single-stage problems to

470

multi-stage and multiple-shooting OCP (1) is straightforward. The adaptive

471

wavelet-based control grid refinement is performed for each stage k ∈ K until

472

the detected switching structure Ak between two successive refinement steps

473

remains unchanged for all stages. In the second step, we formulate a multi-

474

stage problem with underlying algorithmic stages, Eq. (2), which considers

475

the detected switching structure. After refinement, each sensitivity-seeking

476

and path-constrained arc for all physical stages k ∈ K is again analyzed for

478

479

480

481

482

M

d

te

Ac ce p

477

an

469

any structural changes. These structural changes may lead to the introduction or elimination of algorithmic stages l ∈ Lk , ∀k ∈ K. Generally, the solution of the discrete problem involves the discretization

of a continuous OCP by either shooting approaches or full discretization. In application, the user is interested in solving a discrete problem that is a close approximation of the continuous one. Hence, it is desirable that the stopping

483

criterion reflects this property and simultaneously terminates the adaptive

484

refinement procedure after a finite number of refinement iterations.

485

The stopping criterion (10) is easy to evaluate but may lead to better

486

than optimal objective values for path-constrained OCPs. This behavior 27

Page 27 of 72

stems from converting the originally continuous state path-constraint of the

488

OCP into a discrete set of constraints (cf. Eq. (6)). Hence, the information

489

about the continuous nature of the state path-constraint is unknown to the

490

employed NLP solver. Consequently, the solver considers a solution to be su-

491

perior if an improvement in the objective function is possible without violat-

492

ing the discretized state path-constraint. Often this improvement is achieved

493

by introducing constraint violations in the intervals between the grid points

494

of the state path-constraint discretization (cf. Eq. (6)), which are not covered

495

by the NLP problem. This behavior is easily observed by simulating the un-

496

derlying DAE system using some optimal decision variables plk and simulta-

497

neously evaluating the state path-constraints on a grid with a finer resolution

498

than the one employed in the NLP formulation. An exemplary state trajec-

499

tory and state path-constraint showing the described phenomena is given in

500

Figure 2. Hence, it is desirable to minimize these intermediate constraint

501

violations (ICV) and simultaneously improve the objective function. How-

502

ever, a stopping criterion based on the PMP is most rigorous, but requires

508

Ac ce p

te

d

M

an

us

cr

ip t

487

509

of the control variable is finer than that of the state path-constraints or not

510

assessable if the discretization of the control variables is decoupled from the

511

discretization of the state path-constraints.

503

504

505

506

507

the (expensive) solution of the adjoint system (Hannemann and Marquardt, 2012).

Here, we combine the stopping criterion on the objective function (10)

together with checking ICVs. Typically, these ICVs are large for coarsely discretized state path-constraints, even if the optimal solution of the NLP is solved to a high accuracy. This eﬀect may be amplified if the discretization

28

Page 28 of 72

1.2 1.1

ip t

1 state

0.9 0.8 0.7

cr

0.6 0.5 0.2

0.4 0.6 time [s]

0.8

1

us

0.4 0

Figure 2: Left: Typical state trajectory showing constraint violations at intermediate time intervals (hatched surface). Right: An exemplary discretization of a state path-constraint

an

and the associated subgrid points required for the evaluation of the ICV denoted with ◦. In both figures, × markers denote the discretization grid of the state path-constraint,

M

while the dashed line indicates the upper bound.

The ICV can be eﬀectively measured by approximating the surface be-

513

tk−1 tween the state trajectory and the constraint (hatchedtk area in Figure 2).

514

Hence, the state path-constraint is evaluated on an extended set ∆Nkl ⊂

515

∆hlk ⊂ [tlk−1 , tlk ] of grid points, which covers the interval between the grid

516

points of the discretized state path-constraints hk (·). Here, Nh additional

518

519

520

te

Ac ce p

517

d

512

subgrid points are introduced equidistantly in the interval (tlk,j , tlk,j+1 ) with j ∈ ∆hlk and ∀k ∈ K (cf. Figure 2 for an exemplary discretization). Subsequently, the DAE system is evaluated on the refined set ∆Nkl as

follows:

) ( ˜ lk (τ ; plk ) , Mk x˙ lk = fk xlk (τ ; plk ), zlk (τ ; plk ), u ) ( ˜ lk (τ ; plk ) , 0 = gk xlk (τ ; plk ), zlk (τ ; plk ), u ] [ τ ∈ tlk−1 , tlk ∪ ∆Nkl , ∀k ∈ K.

(11a) (11b) (11c)

29

Page 29 of 72

This allows to estimate the violation of an active constraint hk (·)(cf. Figure

522

2, hatched surface). Hence, using the Riemann sum the ICV of constraint

523

hk (·) can be approximated as

ip t

521

an

us

cr

∑

hk (tk,j ) (tk,j+1 − tk,j ) ω (hk (tk,j )) , ϵk,I =

j∈∆Nkl 1   1 for hk (tk,j ) < 0, ω(hk (tk,j )) = ,  0 else.

(12b)

M

∀k ∈ K.

(12a)

This measure allows judging the quality of an optimal solution by a simple

525

forward simulation of the DAE system. Further, it can be applied as an ad-

526

ditional stopping criterion for an adaptive refinement approach by requiring

527

that

529

530

531

te

Ac ce p

528

d

524

∑

ϵk,I ≤ ϵI ,

(13)

k

where ϵI denotes a user-specified feasibility tolerance. However, the com-

putational overhead required for its evaluation depends on the size of the underlying process model but it may be cheaply evaluated in many cases. 4. Implementation and solvers

532

The presented algorithms have been implemented in DyOS (Dynamic

533

Optimization Software), a software package developed and maintained by

534

AVT - Process Systems Engineering, RWTH Aachen University.

30

Page 30 of 72

Algorithm 1: Adaptive multi-stage switching structure detection.

ip t

(Continued in Algorithm 2) choose maximum number of refinement steps: λmax , κmax ,

2

choose stopping criterion parameters (cf. Eqs. (10) and (12)):

cr

1

Nkl , ϵI , ϵΦ

4

if multiple-shooting problem then

us

3

choose number of physical stages / multiple-shooting intervals K end if

6

formulate multi-stage OCP according to Eq. (1)

7

choose approximation order r and discretization grid ∆ulk,i for control

an

5

M

variables u˜lk,i in Eq. (3)

discretize state path-constraints in Eq. (6)

9

set λ = 0

d

8

initialize decision variables pl,λ k in Eq. (7) or Eq. (8)

11

set algorithmic stages Lk = 1 and Aλk = {sens} control grid adaptation loop.

Ac ce p

/* Step 1:

te

10

12

*/

while (λ < λmax ) do

13

solve NLP problem in Eq. (7) or Eq. (8) for pl,λ k

14

analyze solution structure for all stages k ∈ K → Aλ+1 (cf.

15

Section 3.5) ( ) if Aλ ̸= Aλ+1 & (λ < λmax ) then

16

apply wavelet-based control grid refinement (cf. Section 3.4)

17

discretize state path-constraints in Eq. (6)

18

update vector pl,λ+1 in Eq. (7) or Eq. (8) using current k decision variables pl,λ k

19 20

else go to line 17

21

end if

22

set λ = λ + 1

23

end while

31

Page 31 of 72

Algorithm 2: Adaptive multi-stage switching structure detection.

25

set κ = 0

26

set Aκ = Aλ

27

for κ = 0, 1, . . . , κmax do

28

set up Lk algorithmic stages for multi-stage problem (2) using

an

switching structure Aκ

solve multi-stage NLP problem in Eq. (7) or Eq. (8) for pl,κ k

30

analyze solution structure for all stages k ∈ K → Aκ+1 (cf. Section 3.5) if (Aκ ≡ Aκ+1 ) then

M

29

31

*/

cr

structure detection loop with adaptation.

us

/* Step 2:

ip t

(Continuation from Algorithm 1)

evaluate objective improvement in Eq. (10)

33

evaluate ICV in Eq. (12)

34

if (objective improvement > ϵΦ ) & (constraint violation ≤ ϵI )

te

d

32

Ac ce p

& (for κ > 0) then exit

35

36

/* optimal solution found.

*/

else

apply wavelet-based control grid refinement (cf. Section

37

3.4)

38

discretize state path-constraints in Eq. (6)

39

update vector pl,κ+1 in Eq. (7) or Eq. (8) using current k decision variables pl,κ k

40

end if

41

end if

42

set κ = κ + 1

43

end for

32

Page 32 of 72

Note that the objective, constraint and gradient evaluation required for

536

optimization is obtained from simultaneous numerical integration of the

537

DAE system (Eqs. (1b)-(1c)) together with the first-order sensitivity sys-

538

tem (Schlegel et al., 2004; Vassiliadis et al., 1994a,b). The reformulated NLP

539

summarized in Eqs. (7) and (8) are solved by sequential quadratic program-

540

ming (SQP) (Nocedal and Wright, 1999). SNOPT employs a sparse SQP

541

algorithm that does not specifically exploit the block structure resulting

542

from the discretization using multiple-shooting. However, we believe that

543

the use of an oﬀ-the-shelf sparse SQP leverages most of the potential regard-

544

ing computational eﬃciency to solve the NLP resulting from the application

545

of multiple-shooting. Tailored NLP solvers (Bock and Plitt, 1984), which

546

exploit the problem-specific block structure, are expected to oﬀer some ad-

547

ditional computational benefit, but their implementation in the context of

548

adaptive discretization is involved.

te

d

M

an

us

cr

ip t

535

The illustrative examples considered in Section 5 are solved using the NLP

550

solver SNOPT (Gill et al., 1998) and S-LIMEX (Schlegel et al., 2004), which

551

552

553

554

555

Ac ce p

549

is an extension of LIMEX (Deuflhard et al., 1987; Deuflhard and Nowak, 1987), for numerical integration. The models of both examples have been implemented in the modeling environment Modelica8 . The model derivatives are provided using the JADE framework (Hannemann et al., 2010), which applies automatic diﬀerentiation. 8

http://www.modelica.org/

33

Page 33 of 72

556

5. Illustrative example We apply the multi-stage automatic structure detection to the Williams-

558

Otto (WO) semi-batch reactor as introduced by Forbes (1994) based on

559

the original work of Williams and Otto (1960) and thoroughly analyzed by

560

Hannemann and Marquardt (2012). In the stirred tank reactor considered,

561

three exothermic reactions

an

A + B → C,

us

cr

ip t

557

C + B → P + E,

M

P + C → G,

occur in a homogenous, ideally mixed liquid phase. The heat produced by

563

these reactions is removed from the reactor by a cooling jacket, whose cooling

564

capacity is controlled through the cooling water inlet temperature TW (t).

566

567

568

569

570

571

te

The model equations and parameters are given in Appendix A. The pro-

Ac ce p

565

d

562

cess model consists of nx = 9 diﬀerential and nz = 3 algebraic state as well as nu = 3 control variables, i.e., the cooling water temperature TW (t) and the feed flow rates FA (t) and FB (t) of reactants A and B, respectively. In contrast to Hannemann and Marquardt (2012), we have included an additional control variable, i.e. the feed rate FA (t), to mimic a multi-stage recipe control. We do not intend to formulate a coherent recipe control problem from

572

a process engineering perspective, rather, we chose the formulation such that

573

it covers a substantial part of the capabilities of the presented algorithm.

574

First, we mimic a multi-stage recipe control problem with a chang-

575

ing set of controls and constraints. Second, a multiple-shooting example 34

Page 34 of 72

is introduced using the original model and optimal control formulation by

577

Hannemann and Marquardt (2012).

578

5.1. Multi-stage optimal control problem

ip t

576

A multi-stage optimization problems consisting of K = 2 stages that

580

mimics a multi-stage recipe OCP is introduced. The amount of product P in

581

the batch by is maximized according to

Φ(t2 ) dΦ = −ρP V, dt

(14a) Φ(t0 ) = 0,

M

s.t.

an

min FA (t),FB (t),TW (t),∆t1 ,∆t2

us

cr

579

process model equations (A.1).

(14b) (14c)

The problem is formulated to match Eq. (1). The stage lengths ∆t1 and

583

∆t2 are decision variables of the OCP. Both stages are subject to the process

584

model and the corresponding initial conditions defined in Eqs. (A.1).

586

te

Additionally, stage k = 1 defines control and state path-constraints as

Ac ce p

585

d

582

follows

0 dm3 /s ≤ FA (t) ≤ 5.784 dm3 /s, 20◦ C ≤ TW (t) ≤ 100◦ C,

∀t ∈ [t0 , t1 ],

(15a)

∀t ∈ [t0 , t1 ],

(15b)

10−4 s ≤ ∆t1 ≤ 1000 s,

(15c)

60 ◦ C ≤ TR (t) ≤ 70 ◦ C,

∀t ∈ [t0 , t1 ],

70 ◦ C ≤ TR (t1 ) ≤ 70 ◦ C.

(15d) (15e)

35

Page 35 of 72

Feed rate FB (t) is not a control variable in the first stage, where it is fixed to

588

FB (t) = 0 dm3 /s. The remaining decision variables are constrained by Eqs.

589

(15a) - (15c). Further, a path and an endpoint constraint for the temperature

590

are defined in Eqs. (15d) and (15e), respectively.

cr

ip t

587

The second stage reintroduces the feed rate FB (t) as a decision variable

592

but eliminates the feed rate FA (t) and fixes it to 0 dm3 /s. The control and

593

state path-constraints are defined as

an

us

591

0 dm3 /s ≤ FB (t) ≤ 5.784 dm3 /s,

M

20◦ C ≤ TW (t) ≤ 100◦ C,

∀t ∈ [t1 , t2 ],

(16a)

∀t ∈ [t1 , t2 ],

(16b)

10−4 s ≤ ∆t2 ≤ 1000 s,

∀t ∈ [t1 , t2 ],

d

70 ◦ C ≤ TR (t) ≤ 90 ◦ C,

(16c)

0 m3 ≤ V (t2 ) ≤ 5 m3 .

595

596

597

598

599

te

(16e)

The previously introduced endpoint constraint Eq. (15e) for reactor temper-

Ac ce p

594

(16d)

ature TR (t) matches the lower bound of the path-constraint for the reactor temperature of the second stage Eq. (16d). Here, the decision variables are constrained by Eqs. (16a) - (16c). An endpoint constraint for the reactor volume V (t) is specified in Eq. (16e). Finally, the time horizon of the OCP is constrained as

200 s ≤ ∆t1 + ∆t2 ≤ 2000 s, ∆tk = tk+1 − tk ,

(17a) ∀k ∈ K.

(17b)

36

Page 36 of 72

Table 2: Computational statistics for illustrative multi-stage OCP.

λ

constr. vio.

dec. vars.

maj. iter.

time [s]

Step 1: control grid adaptation loop. -1039.27

1.0899

70

115

2

-1039.85

0.5008

70

55

3

-1040.02

0.0976

84

43

κ

time [s]

15.66

15.66

8.00

23.66

7.74

31.40

us

1

∑

ip t

obj. function

cr

ref. iter.

Step 2: adaptive structure detection loop. -1040.12

1.1670

23

82

5.46

36.86

2

-1040.12

0.1493

31

35

2.71

39.57

3

-1040.12

0.0172

41

7

0.70

40.27

4

-1040.12

0.0028

63

3

0.44

40.71

M

an

1

Finely discretized equidistant control grid.

601

602

603

604

605

606

d

0.0384

262

80

49.41

49.41

5.2. Numerical results for the multi-stage optimal control problem The various settings for the algorithmic parameters and tolerances em-

Ac ce p

600

-1040.05

te

eq.

ployed in this example are summarized in Tables B.1 and B.2. The multi-stage problem is initially discretized by introducing a single

algorithmic stage for each physical stage k ∈ K resulting in L1 = L2 = 1. This initialization strategy is reasonable for problems with unknown control structure. Furthermore, we discretize each algorithmic stage l ∈ L using

607

16 equidistant piecewise linear basis functions. This dyadic grid allows the

608

wavelet-based adaptation to start on a low wavelet-scale and avoids deletion

609

of grid points at an early stage of the algorithm.

610

The state path-constraints introduced in Eqs. (15d) - (16d) are discretized 37

Page 37 of 72

according to Eq. (6) (cf. Section 3.1). Thus, the constraint grids consist of

612

32 grid points for both stages. Further, we initialize the duration of the first

613

stage as ∆τ11 = ∆t1 = 100.0 s and of the second stage as ∆τ21 = ∆t2 = 900.0

614

s. This results in 38 decision variables for the first NLP. Table B.4 summarizes

615

the initial values of the decision variables used for the first optimization run.

616

The control grid adaptation loop (Step 1) terminates to a non-changing

617

control switching structure after λ = 3 refinement iterations (see Table C.1).

618

The subsequent structure detection loop (Step 2) fulfills the stopping crite-

619

ria Eq. (10) and Eq. (12) after κ = 4 refinement iterations, while having

620

the same control switching structure in two subsequent iterations. Table 2

621

summarizes the computational statistics for both the control grid adapta-

622

tion loop (Step 1) and the structure detection loop (Step 2), respectively.

623

Further, the table gives a finely discretized equidistant solution for compar-

624

ison. Its columns show the number of refinement iterations, the objective

625

function, the ICV, the number of decision variables, the number of major

626

iterations required by the SQP solver, the time required per iteration, and

628

629

630

631

632

cr

us

an

M

d

te

Ac ce p

627

ip t

611

the accumulated time.

The evolution of the control profiles for the control grid adaptation (Step

1) and the structure detection loop (Step 2) are shown in Figures 3 and 5, respectively. Figures 4 and 6 show exemplarily the path-constraint violation for physical stage k = 1 at refinement iterations λ = 1, and λ = 3 as well as κ = 1, and κ = 4, respectively.

633

We observe for all refinement iterations λ and κ that the total time con-

634

straint formulated in Eq. (17a) is active at its upper bound. Hence, the stage

635

duration constraints (15c) and (16c) are also active.

38

Page 38 of 72

636

5.2.1. Control grid adaptation loop The development of the identified switching structure is depicted in Table

638

C.1. The algorithm identifies a rich switching structure for feed flow rate

639

FA (t) consisting of 5 arcs already in the first refinement iteration λ = 1. In

640

the second refinement step, the increased control grid resolution leads to the

641

evolution of a path-constrained arc with an active lower bound (cf. Figure

642

4). If only the feed flow rate FA (t) was regarded, one may mistake this arc

643

as numerical oscillation, since the control value changes aggressively on a

644

small time horizon. The oscillatory behavior observed for the feed flow rate

645

FA (t) approximately at time point 500 s (cf. Figure 3, bottom left) stems

646

from the deficiency of the adaptive control grid adaptation to exactly catch

647

the switching time between the two arcs. Such behavior has been observed

648

before for the wavelet-based refinement approach (Schlegel et al., 2005). The

649

final adaptation step λ = 3 does not add any novel structural information

650

but increases the confidence that the correct structure has been identified.

cr

us

an

M

d

te

The optimal solution of feed rate FB (t) consists of two intervals (cf. Table

657

Ac ce p

651

ip t

637

658

zon (cf. Table C.1). The initially equidistant parameterization is not suitable

659

to exactly determine the switching structure. The structural information re-

660

mains unchanged from refinement step λ = 2 to λ = 3. The control grid

652

653

654

655

656

C.1). The first arc is determined by an active path-constraint in the first part of the second physical stage. Finally, the control is at its lower bound until the end of the time horizon. This structure remains unchanged over the course of adaptation steps λ. The cooling water temperature TW (t) consists of an altering bang-bang

structure, which ends in a sensitivity-seeking arc at the end of the time hori-

39

Page 39 of 72

2 1

3 2 1

0

0

−1 0

−1 0

500

1000 time [s]

1500

2000

ip t 80

us

3

4

500

60 40 20

an

3

4

100

cr

5

W

6

5

cooling water temperature T [°C]

6

feed rate FB [m /s]

3

feed rate FA [m /s]

adaptation step λ = 1

1000 time [s]

1500

2000

0

500

1000 time [s]

1500

2000

500

1000 time [s]

1500

2000

500

1000 time [s]

1500

2000

3

4 3 2

2 1 0

500

1000 time [s]

−1 0

te

0 −1 0

3

d

1

4

1500

2000

500

1000 time [s]

1500

100

W

5

cooling water temperature T [°C]

6

5

M

6

feed rate FB [m /s]

3

feed rate FA [m /s]

adaptation step λ = 2

2000

80 60 40 20 0

5

3

4 3 2 1

4 3 2 1

0

0

−1 0

−1 0

500

1000 time [s]

1500

2000

100

W

6

5

cooling water temperature T [°C]

6

feed rate FB [m /s]

3

feed rate FA [m /s]

Ac ce p

adaptation step λ = 3

500

1000 time [s]

1500

2000

80 60 40 20 0

Figure 3: Control profiles for adaptation steps λ =1,2 and 3, respectively. The vertical dashed line represents the boundary between the physical stages.

40

Page 40 of 72

66

69.9

64

60.1

60 58 0

200

400 600 time [s]

800

1000

60 59.9

an

R

°

T [ C]

62

600

700 time [s]

800

cr

°

R

70

us

T [ C]

68

R

temperature T [°C]

70

ip t

70.1

72

120

140 time [s]

160

M

Figure 4: Evaluation of Eq. (11) for the temperature constraint for stage k = 1 and adaptation refinement steps λ = 1 (♢), and λ = 3 (·), respectively. The top-right and bottom-right figures represent cut-outs that clarify the ICVs at the upper and lower bound,

d

respectively.

refinement of the cooling water temperature TW (t) improves the recognition

662

of the bang-bang structure.

664

665

666

667

668

Ac ce p

663

te

661

Generally, the control grid refinement loop (Step 1) requires less com-

putational eﬀort compared to a finely resolved equidistant control grid (Schlegel et al., 2005) and significantly improves the recognition of the control switching structure. Here, the number of major iterations required by the NLP solver as well as the time per refinement iteration (cf. Table 2, Step 1, columns 5 and 6) drop significantly from refinement iteration λ = 1 to

669

λ = 2. This results from a good initialization of the decision variables in

670

refinement step λ = 2 using the optimal solution from refinement step λ = 1.

671

Refinement iteration λ = 3 requires only few major NLP iterations less than

41

Page 41 of 72

λ = 2. However, the time required per iteration remains similar despite the

673

increase in the number of decision variables.

ip t

672

The objective function only improves little during the control grid adap-

675

tation loop (cf. Table 2, Step 1, column 2), which is an expected behavior

676

(Schlegel et al., 2005). However, the ICV (cf. Table 2, Step 1, column 3)

677

is successfully reduced by the proposed coupling of the path-constraint and

678

the control discretization (cf. Section 3.1). A visual inspection of the reactor

679

temperature path-constraint evaluation (cf. Figure 4, right) confirms that

680

the ICV is reduced by the adaptive control grid adaptation. However, we

681

observe still some overshoot, especially for the active upper bound, due to

682

the deficiency to exactly determine the switching points of the active path-

683

constraints.

684

5.2.2. Adaptive structure detection loop

d

M

an

us

cr

674

The adaptive structure detection loop (Step 2) exploits the previously de-

686

tected switching structure by reparameterizing the multi-stage problem using

688

689

690

691

692

Ac ce p

687

te

685

algorithmic stages. We choose 2 piecewise linear approximations for control arcs with inactive bounds as the initial discretization. This parameterization is refined in the course of the structure detection loop κ. The switching structure of feed rates FA (t), FB (t), and cooling water tem-

perature TW (t) remains unchanged in the course of the structure detection loop (Step 2) (cf. Figure 5). The proposed reformulation does not only im-

693

prove the solution quality, rather it significantly drops the number of decision

694

variables in the first refinement iteration κ = 1 (cf. Table 2, Step 2, column

695

4). The number of major iterations required by the SQP solver increases

696

due to nonlinearities arising from the optimization of the algorithmic stage 42

Page 42 of 72

lengths ∆τkl (cf. Table 2, Step 2, column 5). However, the computational

698

time required for the first iteration κ = 1 is smaller than that for the final

699

control grid adaptation λ = 3. The number of major iterations required by

700

the SQP solver as well as the computational time per refinement iteration

701

κ drops significantly in the course of the structure detection loop (Step 2).

702

The initialization of the decision variables with the previous optimal solution

703

leads to fast convergence.

us

cr

ip t

697

The numerical optimization algorithm is capable of further improving the

705

objective function by exactly catching the switching structure (cf. Table 2,

706

Step 2, column 2). Initially, the ICV increases due to the coarse discretization

707

of the controls at path-constrained arcs (cf. Table 2, Step 2, column 3). But

708

with increasing resolution these constraint violations are further decreased.

709

The previously observed overshooting of the path-constraint is successfully

710

compensated with the proposed structure detection approach (cf. Figure 6,

711

right). Thus, we obtain a clear solution structure in the controls and the

712

path-constraints. We expect that the ICVs disappear with increasing con-

714

715

716

717

M

d

te

Ac ce p

713

an

704

trol grid resolution, if the correct switching structure is detected. Hence, the optimal solution is not only superior in the objective but also in satisfying path-constraints compared to the grid refinement approach and an equidistant discretization. Additionally, the equidistant control grid discretization is inferior in the required computational load.

43

Page 43 of 72

3 2 1 0

3 2 1 0

500

1000 time [s]

1500

−1 0

2000

500

1000 time [s]

1500

80 60

2000

ip t

40 20

0

500

1000 time [s]

1500

2000

500

1000 time [s]

1500

2000

500

1000 time [s]

1500

2000

an

−1 0

4

us

3

4

100

cr

5

W

6

5

cooling water temperature T [°C]

6

feed rate FB [m /s]

3

feed rate FA [m /s]

structure detection step κ = 1

3 2 1 0

3 2 1 0

500

1000 time [s]

1500

2000

−1 0

500

te

−1 0

4

d

3

4

100

W

5

cooling water temperature T [°C]

6

5

M

6

feed rate FB [m /s]

3

feed rate FA [m /s]

structure detection step κ = 2

1000 time [s]

1500

2000

80 60 40 20 0

structure detection step κ = 4

feed rate FB [m /s]

4

3

3

feed rate FA [m /s]

3 2 1

6 5 4 3 2 1

0

0

−1 0

−1 0

500

1000 time [s]

1500

2000

100

W

5

cooling water temperature T [°C]

Ac ce p

6

500

1000 time [s]

1500

2000

80 60 40 20 0

Figure 5: Control profiles for structure detection steps κ = 1,2 and 4, respectively. The vertical dashed line represents the boundary between the physical stage. The vertical

dash-dotted lines indicate the beginning and the end of the underlying algorithmic stages.

44

Page 44 of 72

70 69.9

64

60.1

60 58 0

200

400 600 time [s]

800

1000

60 59.9

an

R

°

T [ C]

62

600

700 time [s]

800

cr

66

us

°

R

68

R

temperature T [°C]

70

ip t

70.1 T [ C]

72

120

140 time [s]

160

M

Figure 6: Evaluation of Eq. (11) for the temperature constraint for stage k = 1 and switching structure refinement steps κ = 1 (♢), and κ = 4 (·), respectively. The top-right and bottom-right figures represent cut-outs that clarify the ICVs at the upper and lower

718

5.3. Multiple-shooting optimal control problem Here, the OCP originally formulated by Hannemann and Marquardt

725

Ac ce p

719

te

d

bound, respectively.

726

the number of shooting intervals has no eﬀect on the automatic switching

727

structure detection. The stage length of each individual stage/shooting in-

728

terval is fixed to ∆tk = 250 s. Operation aims at maximizing the conversion

720

721

722

723

724

(2012) is solved. Assassa and Marquardt (2014) demonstrated the potential of adaptive grid refinement for the WO semi-batch reactor using adaptive multiple-shooting. The same example is used here to demonstrate the possible improvement by adaptive multiple-shooting with structure detection. Here, we arbitrarily introduce K = 4 independent physical stages or

multiple-shooting intervals on the fixed time horizon of 1000 s. Note that

45

Page 45 of 72

of the desired products P and E; it is formulated as a minimization problem

730

as follows

Φ(t4 )

cr

min

FB (t),TW (t),x0k

ip t

729

Φ(t0 ) = 0,

the

originally

introduced

problem

(18b) (18c) (18d)

731

To

formulation

732

(Hannemann and Marquardt, 2012), we have to set the feed rate FA (t) = 0

733

m3 /s for all shooting intervals k ∈ K and exclude it from the optimization

734

problem. All stages k ∈ K are subject to the process model defined in Eqs.

735

(A.1) as well as control and state path-constraints summarized as

te

d

M

meet

an

process model equations (A.1).

us

ME r2 ρB ρC V, s.t. Φ˙ = −cp (k2 r2 ρB ρC − k4 r3 ρC ρP ) V − ce MB

(18a)

Ac ce p

0 = xk (tk−1 ) − xk,0 ,

Up xLo k,0 ≤ xk,0 ≤ xk,0 ,

0 = xk (tk ) − xk+1 (tk ),

0 dm3 /s ≤ FB (t) ≤ 5.784 dm3 /s

∀k ∈ K,

(19a)

∀k ∈ K,

(19b)

∀k ∈ K\K,

(19c)

∀t ∈ [tk−1 , tk ],∀k ∈ K,

(19d)

20◦ C ≤ TW (t) ≤ 100◦ C

∀t ∈ [tk−1 , tk ],∀k ∈ K,

(19e)

60 ◦ C ≤ TR (t) ≤ 90 ◦ C

∀t ∈ [tk−1 , tk ],∀k ∈ K.

(19f)

736

Here, x = [ρA , ρB , ρC , ρE , ρG , ρP , TR , V, Φ]T refers to the diﬀerential states of

737

the process model and the objective function (18b). Additionally, we define

738

an endpoint constraint on the volume for the final stage k = 4 as 46

Page 46 of 72

V (t4 ) ≤ 5 m3 .

ip t

(20)

The multiple-shooting OCP can be summarized by Eqs. (18) - (20).

740

5.4. Numerical results for multiple-shooting optimal control problem

742

us

741

cr

739

The settings for the algorithmic parameters and tolerances employed in this example are summarized in Tables B.1 and B.3.

The multiple-shooting problem introduced in Section 5.3 is initially dis-

744

cretized by introducing a single algorithmic stage for each physical stage

745

resulting in Lk = 1. Furthermore, we discretize each algorithmic stage l ∈ L

746

using 1 equidistant piecewise linear basis function, which corresponds to the

747

simplest possible discretization. The control approximations are discontin-

748

uous at numerical and physical stage boundaries. All control values are

749

initially chosen such that both control variables TW (t) and FB (t) are at their

750

upper bounds (cf. Table B.5), respectively. The path-constraint introduced

752

753

754

755

756

M

d

te

Ac ce p

751

an

743

in Eq. (19f) is discretized on a grid twice as fine as the control grid (cf. Section 3.1).

Table 3 summarizes the computational statistics for both the control grid

adaptation loop (Step 1) and the structure detection loop (Step 2), respectively. Further, the table gives a finely discretized equidistant solution for comparison. The control grid adaptation loop (Step 1) terminates to a non-

757

changing control switching structure after λ = 4 refinement iterations (cf.

758

Table C.2). The subsequent structure detection loop (Step 2) fulfills both

759

stopping criteria (10) and (12) after κ = 5 refinement iterations. Here, the

47

Page 47 of 72

ICV tolerance is set to ϵI = 10−4 , which is smaller compared to the previous

761

example (cf. Table B.3).

ip t

760

The evolution of the control profiles for the control grid adaptation loop

763

(Step 1) and the structure detection loop (Step 2) are shown in Figures 7 and

764

8, respectively. The path-constraint plot was omitted since a similar behavior

765

in the decrease of ICV was already observed in the previous example. Note

766

that the quasi-analytical solution of the WO semi-batch reactor has been

767

previously published (Hannemann and Marquardt, 2012). The control and

768

state trajectories obtained with the presented approach are nearly identical

769

compared to the quasi-analytical solution.

770

5.4.1. Control grid adaptation loop

M

an

us

cr

762

The feed flow rate FB (t) is already reasonably well approximated in re-

772

finement iteration λ = 1 (cf. Figure 7, left). The control switching structure

773

shows a bang-bang behavior switching from the upper to the lower bound at

774

approximately half of the time horizon. The adaptive control grid refinement

776

777

778

779

780

te

Ac ce p

775

d

771

is not able to exactly resolve the switching point even with increasing refinement steps λ. This is typical for the wavelet-based refinement approach (Schlegel et al., 2005) and has also been observed for this problem before (Assassa and Marquardt, 2014). At refinement step λ = 1 the cooling water temperature TW (t) consists of

sensitivity-seeking arcs for all physical stages k ∈ K (cf. Figure 7, right). The

781

improved control grid resolution in refinement step λ = 2 provides already

782

a reasonable approximation for the second half of the control trajectory.

783

However, the control structure is not clearly identifiable for the first half of

784

the trajectory. Finally, the control switching structure remains unchanged 48

Page 48 of 72

3 2 1 0 400 600 time [s]

800

1000

60 40 20 0

200

400 600 time [s]

800

1000

400 600 time [s]

800

1000

400 600 time [s]

800

1000

an

200

80

cr

4

us

3

feed rate FB [m /s]

5

−1 0

100

W

cooling water temperature T [°C]

6

ip t

adaptation step λ = 1

M

4 3

d

2 1 0

400 600 time [s]

Ac ce p

200

te

3

feed rate FB [m /s]

5

−1 0

800

1000

40 20 200

100

W

4 3 2 1 0

400 600 time [s]

60

0

cooling water temperature T [°C]

3

feed rate FB [m /s]

5

200

80

adaptation step λ = 4

6

−1 0

100

W

6

cooling water temperature T [°C]

adaptation step λ = 2

800

1000

80 60 40 20 0

200

Figure 7: Control profiles for adaptation steps λ =1, 2, and 4, respectively. The vertical dashed line represents the boundary between the physical stages.

49

Page 49 of 72

Table 3: Computational statistics for illustrative multiple-shooting OCP.

λ

constr. vio.

dec. vars.

maj. iter.

time [s]

Step 1: control grid adaptation loop. -4756.40

37.177

43

77

2

-4766.71

0.5225

64

38

3

-4767.74

0.0734

73

34

4

-4768.09

0.0130

85

11.01

2.53

13.54

2.64

16.18

1.36

17.54

Step 2: adaptive structure detection loop.

an

κ

15

time [s]

11.01

us

1

∑

ip t

obj. function

cr

ref. iter.

-4768.31

0.0656

43

43

2.28

19.82

2

-4768.31

0.0081

47

9

0.92

20.74

3

-4768.31

0.0007

55

3

0.39

21.13

4

-4768.31

0.0002

71

4

0.42

21.55

5

-4768.31

0.0000

95

1

0.14

21.69

d

M

1

785

786

787

788

789

790

-4768.09

0.0056

163

82

33.07

33.07

Ac ce p

eq.

te

Finely discretized equidistant control grid.

and is identified correctly at refinement step λ = 4. The control grid refinement loop (Step 1) requires only little com-

putational eﬀort compared to a finely resolved equidistant control grid (Schlegel et al., 2005) and significantly improves the recognition of the control switching structure as observed for the previous example and previously for this example by Assassa and Marquardt (2014). Here, the number of ma-

791

jor NLP iterations but especially the time per refinement iteration (cf. Table

792

3, Step 1, columns 5 and 6) drop significantly from refinement iteration λ = 1

793

to λ = 2. Refinement iterations λ = 3 and 4 require only few major NLP

50

Page 50 of 72

iterations less than refinement iteration λ = 2. The time required per refine-

795

ment iteration remains similar, despite the increase in the number of decision

796

variables.

ip t

794

The objective function only improves little during the control grid adap-

798

tation loop (cf. Table 3, Step 1, column 2) which is an expected behavior

799

already observed in the previous example. However, the ICV (cf. Table 3,

800

Step 1, column 3) is successfully reduced by the proposed coupling of the

801

path-constraint and the control discretization (cf. Section 3.1).

802

5.4.2. Adaptive structure detection loop

an

us

cr

797

The adaptive structure detection loop (Step 2) exploits the previously de-

804

tected switching structure by reparameterizing the multiple-shooting problem

805

using algorithmic stages. Here, we choose 2 piecewise linear approximations

806

for control arcs with inactive bounds as the initial discretization. This pa-

807

rameterization is refined in the course of the structure detection loop κ.

809

810

811

812

813

814

d

te

The switching structure of feed rate FB (t), and cooling water temperature

Ac ce p

808

M

803

TW (t) remains unchanged in the course of the structure detection loop (Step 2) (cf. Figure 8). The proposed reformulation leads in the first refinement step κ = 1 to the theoretical optimal objective function of Φ = −4768.31 (Hannemann and Marquardt, 2012) with only a moderate increase in the ICV (cf. Table 3, Step 2, column 2). The computational load is still comparable to the final refinement step of the control grid adaptation λ = 4 although

815

the number of major NLP iterations increased (cf. Table 3, Step 2, column

816

5). The number of major NLP iterations as well as the computational time

817

per refinement iteration κ drops significantly in the course of the structure

818

detection loop (Step 2). 51

Page 51 of 72

structure detection step κ = 1

ip t

3 2 1 0 200

400 600 time [s]

800

1000

80 60

cr

4

40 20 0

us

3

feed rate FB [m /s]

5

−1 0

100

W

cooling water temperature T [°C]

6

200

400 600 time [s]

800

1000

800

1000

800

1000

80

M

4 3 2

d

1 0 200

te

3

feed rate FB [m /s]

5

−1 0

100

W

6

cooling water temperature T [°C]

an

structure detection step κ = 3

400 600 time [s]

800

1000

60 40 20

0

200

400 600 time [s]

Ac ce p

structure detection step κ = 5

3

feed rate FB [m /s]

5 4 3 2 1 0

−1 0

200

400 600 time [s]

100

W

cooling water temperature T [°C]

6

800

1000

80 60 40 20 0

200

400 600 time [s]

Figure 8: Control profiles for structure detection steps κ = 1, 3, and 5, respectively. The vertical dashed line represents the boundary between the physical stage. The vertical dash-dotted lines indicate the beginning and the end of the underlying algorithmic stages.

52

Page 52 of 72

The ICV (cf. Table 3, Step 2, column 3) becomes insignificant, which re-

820

inforces the claim that the heuristic criterion to discretize the path-constraint

821

(cf. Section 3.1) is well chosen. Thus, we obtain a clear and accurate switch-

822

ing structure in the controls and the path-constraints. The optimal solution

823

is superior in the objective, in satisfying the constraints, and computational

824

load compared to the equidistant discretization (cf. Table 3).

825

6. Conclusions

an

us

cr

ip t

819

This paper generalizes the adaptive switching structure approach intro-

827

duced by Schlegel and Marquardt (2006b) from single-stage problems and

828

single-shooting to multi-stage problems and multiple-shooting. The appli-

829

cation of multi-stage OCP allows formulating so called physical stages with

830

a varying set of controls, objectives, constraints, and model equations, re-

831

spectively. The proposed formulation allows to regard the physical stages

832

as multiple-shooting intervals by including the initial values of all physical

833

stages to the decision-variable vector.

835

836

837

838

839

d

te

Ac ce p

834

M

826

The extended approach is divided in mainly two steps: The first step gen-

erates a problem-dependent discretization using the wavelet-based refinement approach (Schlegel et al., 2005; Assassa and Marquardt, 2014). An unnecessarily fine control grid is avoided by terminating the refinement procedure if the control switching structure remains unchanged in two consecutive iterations. The second step reformulates the multi-stage problem incorporating

840

the detected switching structure using algorithmic stages. The wavelet-based

841

grid refinement is applied to the nested multi-stage problem until the stop-

842

ping criteria on the objective and the ICVs are fulfilled. A change in the 53

Page 53 of 72

switching structure leads to a reformulation of the underlying problem with

844

subsequent control grid refinement.

ip t

843

This approach is demonstrated for a multi-stage and multiple-shooting

846

OCP using the WO semi-batch reactor. The proposed algorithm successfully

847

detects the control switching structure in both cases. A solution of user-

848

specified accuracy can be obtained using only few decision variables, while

849

ICVs are decreased.

us

cr

845

Further, the optimal control switching structure increases process in-

851

sight rather than just improving the objective function numerically. For

852

example, the insight gained can be used to apply NCO tracking control

853

(Srinivasan et al., 2003a) for non-linear processes by using the solution of

854

an oﬀ-line optimization only (Kadam et al., 2007).

M

an

850

The adaptive structure detection approach for multi-stage and multiple-

856

shooting problems, though powerful in concept, oﬀers still room for improve-

857

ment to even further improve its robustness and eﬃciency.

te

Future work should therefore address the following issues: The wavelet-

864

Ac ce p

858

d

855

865

tackled by choosing suitable scalings for the control variables.

859

860

861

862

863

based control grid refinement may fail in those cases where the control variables show only small changes in their value. In such cases, not yet observed in the cases studies of this paper, it is likely that the control grid is not parameterized eﬃciently. Furthermore, the adaptive refinement procedure assumes that a large frequency change has a large impact on the objective function, which is not necessarily true. Both problems could be partially

866

A diﬀerent approach that addresses the control grid refinement requires

867

substituting the purely signal-based adaptation strategy by a more rigorous

54

Page 54 of 72

strategy which is based on a sensitivity analysis of the parameterized OCP

869

(see Binder et al. (2000) for a first attempt). Composite adjoints can be

870

interpreted as derivatives of the objective function with respect to the deci-

871

sion variables Hannemann and Marquardt (2010), thus, hey are suitable for

872

a sensitivity-based grid refinement.

cr

ip t

868

At this point, the stopping criterion in the structure detection algorithm

874

is of an empirical nature, because it only compares the improvement in the

875

objective function in two subsequent iterations as well as checking the ICV.

876

Further, we cannot guarantee that the ICVs decrease over the course of

877

the switching structure refinement although the ICV was successfully de-

878

creased in the presented examples. In case of constrained optimization, a

879

higher resolution of the controls might improve the objective function but

880

introduce larger ICVs. This undesirable eﬀect could be remedied in future

881

work, if a more rigorous stopping criterion based on the PMP were employed.

882

Hannemann and Marquardt (2012) construct a measure of the distance be-

883

tween the true solution (satisfying PMP) and the approximate solution of the

885

886

887

888

889

890

an

M

d

te

Ac ce p

884

us

873

NLP. Combining this measure with adaptive structure detection algorithm should result in a justifiable indicator to terminate the algorithm. 7. Acknowledgments

This work has been partially funded by the German Research School for

Simulation Sciences.

F. Assassa and W. Marquardt. Dynamic optimization using adaptive direct multiple shooting. Comput. Chem. Eng., 60:242–259, 2014.

55

Page 55 of 72

M. Avraam, N. Shah, and C. Pantelides. Modelling and optimisation of

892

general hybrid systems in the continuous time domain. In European Sym-

893

posium on Computer Aided Process Engineering (ESCAPE 8), pages 221–

894

228. Pergamon-Elsevier Science Ltd., 1998.

cr

ip t

891

R. Bachmann, L. Br¨ ull, T. Mrziglod, and U. Pallaske. On methods for re-

896

ducing the index of diﬀerential-algebraic equations. Comput. Chem. Eng.,

897

14:1271 – 1273, 1990.

us

895

E. Balsa-Canto, J. R. Banga, A. A. Alonso, and V. S. Vassiliadis. Eﬃcient

899

Optimal Control of Bioprocesses Using Second-Order Information. Indus-

900

trial & Engineering Chemistry Research, 39(11):4287–4295, 2000.

903

M

processes. AIChE J., 40(6):966–979, 1994. L. Biegler.

d

902

P. I. Barton and C. C. Pantelides. Modeling of combined discrete/continuous

Solution of dynamic optimization problems by successive

te

901

an

898

quadratic programming and orthogonal collocation. Comp. Chem. Eng.,

905

8:243–248, 1984.

906

907

908

909

910

Ac ce p

904

L. Biegler. Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes. MOS-SIAM Series on Optimization. Society for Industrial and Applied Mathematics, 2010.

L. T. Biegler and I. E. Grossmann. Retrospective on optimization. Computers & Chemical Engineering, 28(8):1169 – 1192, 2004. ISSN 0098-1354.

911

T. Binder, A. Cruse, C. Villar, and W. Marquardt. Dynamic optimization

912

using a wavelet based adaptive control vector parameterization strategy.

913

Comput. Chem. Eng., 24:1201–1207, 2000. 56

Page 56 of 72

T. Binder, L. Blank, H. Bock, R. Bulirsch, W. Dahmen, M. Diehl, T. Kro-

915

nseder, W. Marquardt, J. Schl¨oder, and O. v. Stryk. Introduction to Model

916

Based Optimization of Chemical Processes on Moving Horizons. In On-

917

line Optimization of Large Scale Systems, pages 295–339. Springer Berlin

918

Heidelberg, 2001a.

cr

ip t

914

T. Binder, L. Blank, W. Dahmen, and W. Marquardt. Iterative algorithms

920

for multiscale state estimation, Part 1: Concepts. J. Optim. Theor. Appl.,

921

111(3):501–527, 2001b.

an

us

919

T. Binder, L. Blank, W. Dahmen, and W. Marquardt. Iterative algorithms

923

for multiscale state estimation, Part 2: Numerical investigations. J. Optim.

924

Theor. Appl., 111(3):529–551, 2001c.

M

922

H. Bock. Randwertproblemmethoden zur Parameteridentifizierung in Syste-

926

men Nichtlinearer Diﬀerentialgleichungen. PhD thesis, University of Bonn,

927

Bonner Mathematische Schriften, 1987.

929

930

931

932

te

Ac ce p

928

d

925

H. Bock and K.-J. Plitt. A multiple shooting algorithm for direct solution of optimal control problems. In Proc. of the 9th IFAC World Congress, Budapest. Pergamon Press, 1984.

D. Bonvin. Optimal operation of batch reactors - a personal view. J. Proc. Contr., 8(5-6):355–368, 1998.

933

S. Br¨ uggemann, J. Oldenburg, P. Zhang, and W. Marquardt. Robust dy-

934

namic simulation of three-phase reactive batch distillation columns. Ind.

935

Eng. Chem. Res., 43(14):3672–3684, 2004.

57

Page 57 of 72

937

A. Bryson and Y.-C. Ho. Applied Optimal Control. Taylor & Francis, Bristol, PA, 1975.

ip t

936

C. B¨ uskens. Optimierungsmethoden und Sensitivit¨ atsanalyse f¨ ur optimale

939

Steuerprozesse mit Steuer- und Zustands-Beschr¨ ankungen. PhD thesis,

940

University of M¨ unster, 1998.

943

944

us

cess systems. AIChE Journal, 33(8):1257–1270, 1987.

an

942

J. Cuthrell and L. Biegler. On the optimization of diﬀerential-algebraic pro-

W. Dahmen. Wavelet and multiscale methods for operator equations. Acta Numerica, pages 55–228, 1997.

M

941

cr

938

C. de Boor. A Practical Guide to Splines. Springer-Verlag, New York, 1978.

946

P. Deuflhard and U. Nowak. Extrapolation integrators for quasilinear im-

947

plicit odes. In P. Deuflhard and B. Enquist, editors, Large Scale Scientific

948

Computing, NATO–ASI Series, pages 37–50. Birkh¨auser, Basel, 1987.

950

951

952

953

te

Ac ce p

949

d

945

P. Deuflhard, E. Hairer, and J. Zugck. One-step and extrapolation methods for diﬀerential-algebraic systems. Numerische Mathematik, 51(5):501–516, 1987.

J. Forbes. Model Structure and Adjustable Parameter Selection for Operations Optimizations. PhD thesis, McMaster University, Hamilton, Canada, 1994.

954

P. Gill, W. Murray, and M. Saunders. SNOPT: An SQP algorithm for

955

large-scale constrained optimization. Technical report, Stanford Univer-

956

sity, Stanford, USA, December 1998.

58

Page 58 of 72

R. Hannemann and W. Marquardt. Continuous and discrete composite ad-

958

joints for the Hessian of the Lagrangian in shooting algorithms for dynamic

959

optimization. SIAM Journal on Scientific Computing, 31(6):4675–4695,

960

2010.

cr

ip t

957

R. Hannemann and W. Marquardt. How to verify optimal controls computed

962

by direct shooting methods? - A tutorial. Journal of Process Control, 22:

963

494–507, 2012.

an

us

961

R. Hannemann, W. Marquardt, U. Naumann, and B. Gendler. Discrete first-

965

and second-order adjoints and automatic diﬀerentiation for the sensitivity

966

analysis of dynamic models . Procedia Computer Science, 1:297–305, 2010.

967

R. Hannemann-Tam´as. Adjoint Sensitivity Analysis for Optimal Control of

968

Non-Smooth Diﬀerential-Algebraic Equations. PhD thesis, RWTH Aachen,

969

Berichte aus der Verfahrenstechnik, 2012.

971

972

973

974

975

d

te

R. Hannemann-Tam´as, D. Munoz, and W. Marquardt. Adjoint sensitiv-

Ac ce p

970

M

964

ity analysis for nonsmooth diﬀerential-algebraic equation systems. SIAM Journal of Scientific Computing (submitted), 2014.

A. Hartwich. Adaptive Numerical Methods for Dynamic Real-time Optimization of Chemical Processes. PhD thesis, RWTH Aachen, Berichte aus der Aachener Verfahrenstechnik, 2013.

976

M. Hazewinkel. Encyclopedia of Mathematics. Springer Netherlands, 2001.

977

ISA. ANSI/ISA-88.01-2010 Batch Control Part 1: Models and terminology,

978

2010. 59

Page 59 of 72

D. Jacobson, M. Lele, and J. Speyer. New necessary conditions of optimality

980

for control problems with state-variable inequality constraints. Journal

981

of Mathematical Analysis and Applications, 35(2):255 – 284, 1971. ISSN

982

0022-247X.

cr

ip t

979

J. Kadam, B. Srinivasan, D. Bonvin, and W. Marquardt. Optimal grade

984

transition in industrial polymerization processes via NCO tracking. AICHE

985

J., 53(3):627 – 639, 2007.

990

991

992

993

994

995

996

997

an

M

989

R. I. Leine and H. Nijmeijer. Dynamics and Bifurcations of Non-Smooth Mechanical Systems. Springer, Berlin, 2004.

d

988

ming problems. Comput. Math. Prog., 15:261–280, 1985.

R. Mattheij and G. Staarink. On optimal shooting intervals. Mathematics

te

987

D. Kraft. On converting optimal control problems into nonlinear program-

of Computation, 42(165):25–40, 1984. S. E. Mattsson and G. S¨oderlind. Index reduction in diﬀerential-algebraic

Ac ce p

986

us

983

equations using dummy derivatives. SIAM J. Sci. Comput., 14(3):677– 692, 1993.

H. Maurer and U. Heidemann. In Optimization and Optimal Control, volume 477 of Lecture Notes in Mathematics, pages 244–260. Springer Berlin Heidelberg, 1975.

998

D. D. Morrison, J. D. Riley, and J. F. Zancanaro. Multiple shooting method

999

for two-point boundary value problems. Commun. ACM, 5(12):613–614,

1000

Dec. 1962. 60

Page 60 of 72

1002

S. Nilchan and C. Pantelides. On the optimisation of periodic adsorption processes. Adsorption, 4(2):113–147, March 1998.

ip t

1001

J. Nocedal and S. Wright. Numerical Optimization. Springer-Verlag, 1999.

1004

C. Pantelides. The consistent initialization of diﬀerential-algebraic systems.

1007

us

1006

SIAM J. Sci. and Stat. Comput., 9(2):213–231, March 1988.

D. Ramkrishna. On modeling of bioreactors for control. Journal of Process Control, 13(7):581 – 589, 2003.

an

1005

cr

1003

S. Sager. Reformulations and algorithms for the optimization of switching

1009

decisions in nonlinear optimal control. Journal of Process Control, 19(8):

1010

1238 – 1247, 2009. ISSN 0959-1524.

M

1008

R. Sargent and G. Sullivan. The development of an eﬃcient optimal control

1012

package. In J. Stoer, editor, Proceedings of the 8th IFIP Conference on

1013

Optimization Techniques, pages 158–168. Springer-Verlag, 1978.

1015

1016

1017

1018

1019

te

Ac ce p

1014

d

1011

M. Schlegel and W. Marquardt. Detection and exploitation of the control switching structure in the solution of dynamic optimization problems. Journal of Process Control, 16(3):275–290, March 2006a.

M. Schlegel and W. Marquardt. Adaptive Switching Structure Detection for the Solution of Dynamic Optimization Problems. Ind. Eng. Chem. Res., 45(24):8083–8094, October 2006b.

1020

M. Schlegel, W. Marquardt, R. Ehrig, and U. Nowak. Sensitivity analysis of

1021

linearly-implicit diﬀerential-algebraic systems by one-step extrapolation.

1022

Appl. Num. Math., 48:83–102, 2004. 61

Page 61 of 72

M. Schlegel, K. Stockmann, T. Binder, and W. Marquardt. Dynamic op-

1024

timization using adaptive control vector parameterization. Computers &

1025

Chemical Engineering, 29(8):1731–1751, 2005.

ip t

1023

A. Schwartz. Theory and implementation of numerical methods based on

1027

Runge-Kutta integration for solving optimal control problems. PhD the-

1028

sis, Dept. of Electrical Engineering and Computer Sciences, University of

1029

California at Berkeley, 1996.

an

us

cr

1026

J. L. Speyer, D. H. Jacobson, D. H. Jacobson, M. M. Lele, and L. Speyerx.

1031

New Necessary Conditions of Optimality for Control Problems with State-

1032

Variable Inequality Constraints. Journal of Mathematical Analysis and

1033

Applications, 35:255–284, 1971.

M

1030

B. Srinivasan, D. Bonvin, E. Visser, and S. Palanki. Dynamic optimization of

1035

batch processes: II. Role of measurements in handling uncertainty. Comp.

1036

Chem. Eng., 27(1):27–44, 2003a.

1038

1039

1040

1041

1042

1043

1044

te

Ac ce p

1037

d

1034

B. Srinivasan, S. Palanki, and D. Bonvin. Dynamic optimization of batch processes I. Characterization of the nominal solution. Comp. Chem. Eng., 27(1):1–26, 2003b.

P. Tanartkit and L. Biegler. A nested, simultanous approach for dynamic optimization problems - II: The outer problem. Comp. Chem. Eng., 21 (12):1365–1388, 1997.

T. Tsang, D. Himmelblau, and T. Edgar. Optimal control via collocation and non-linear programming. Int. J. Control, 21:763–768, 1975.

62

Page 62 of 72

J. Unger, A. Kr¨oner, and W. Marquardt. Structural analysis of diﬀerential-

1046

algebraic equations systems - theory and applications. Comput. Chem.

1047

Eng., 19:867–882, 1995.

ip t

1045

V. Vassiliadis, R. Sargent, and C. Pantelides. Solution of a class of multistage

1049

dynamic optimization problems: 1. Problems without path constraints.

1050

Ind. Eng. Chem. Res., 33(9):2111–2122, 1994a.

us

cr

1048

V. Vassiliadis, R. Sargent, and C. Pantelides. Solution of a class of multistage

1052

dynamic optimization problems: 2. Problems with path constraints. Ind.

1053

Eng. Chem. Res., 33(9):2123–2133, 1994b. O. von Stryk.

M

1054

an

1051

Numerische L¨ osung optimaler Steuerungsprobleme:

Diskretisierung, Parameteroptimierung und Berechnung der adjungierten

1056

Variablen. VDI-Fortschrittsbericht, Reihe 8, Nr. 441. VDI-Verlag, D¨ ussel-

1057

dorf, 1995.

1059

1060

1061

1062

1063

1064

te

T. J. Williams and R. E. Otto. A generalized chemical processing model

Ac ce p

1058

d

1055

for the investigation of computer control. Transactions of the American Institute of Electrical Engineers, Part I: Communication and Electronics, 79(5):458–473, 1960. ISSN 0097-2452.

L. W¨ urth and W. Marquardt. Infinite-Horizon Continuous-Time NMPC via Time Transformation. Automatic Control, IEEE Transactions on, 59(9): 2543–2548, Sept 2014.

63

Page 63 of 72

ip t

The model equations resulting from species and energy balances as well

cr

as reaction kinetics are as follows:

− k1 r1 ρA ρB ,

us

dρA FA ρFA ρA = − (FA + FB ) dt Fscal V Fscal V ρB dρB FB ρFB = − (FA + FB ) dt Fscal V Fscal V − k1 r1 ρA ρB − k2 r2 ρB ρC ,

an

ρA (0) = ρA,0 ,

M

1068

tor

dρC ρC MC =− (FA + FB ) + k1 r1 ρA ρB dt Fscal V MA MC k2 r2 ρB ρC − k3 r3 ρC ρP , − MB dρE ρE ME =− (FA + FB ) + k2 r2 ρB ρC , dt Fscal V MB dρG ρG MG =− (FA + FB ) + k3 r3 ρC ρP , dt Fscal V MC dρP ρP =− (FA + FB ) + k2 r2 ρB ρC dt Fscal V MP − k3 r 3 ρ C ρ P , MC dTR Tin − TR ∆H1 = (FA + FB ) − k1 r1 ρA ρB dt Fscal V CP ∆H3 ∆H2 k2 r2 ρB ρC − k3 r3 ρC ρP − CP CP U Aref − (TR − TW ), CP mref dV FA + FB = , dt Fscal

d

1067

Model equations of Williams-Otto semi-batch reac-

te

1066

Appendix A

Ac ce p

1065

(A.1a)

ρB (0) = ρB,0 ,

(A.1b)

ρC (0) = ρC,0 ,

(A.1c)

ρE (0) = ρE,0 ,

(A.1d)

ρG (0) = ρG,0 ,

(A.1e)

ρP (0) = ρP,0 ,

(A.1f)

TR (0) = TR,0 ,

(A.1g)

V (0) = V0 ,

(A.1h)

64

Page 64 of 72

(

) E1 r1 = exp , TR + 273.15 ( ) E2 r2 = exp , TR + 273.15 ( ) E3 r3 = exp , TR + 273.15

ip t (A.1j)

cr

(A.1k)

us

The parameters and initial values are summarized in Table A.1.

Table A.1: Parameters and initial values of the Williams-Otto semi-batch reactor model.

value unit

reaction rate constant

k2

reaction rate constant

k3

reaction rate constant

∆H1

heat of reaction

∆H2

heat of reaction

∆H3

heat of reaction

E1

activation energy

1.6599E6 m3 /(kg s)−1

te

d

M

k1

an

parameter description

7.2117E8 m3 /(kg s)−1

2.6745E12 m3 /(kg s)−1 -263.8 kJ/kg -158.3 kJ/kg -226.3 kJ/kg -6666.7 K

Ac ce p

1069

(A.1i)

E2

activation energy

-8333.3 K

E3

activation energy

-11111 K

CP

specific heat capacity

cp

revenue for product P

5554.1 $/kg

ce

revenue for product E

125.91 $/kg

Tin

environment temperature

Aref

reference area

9.2903 m2

mref

reference capacity

2105.2 kg

4.184 kJ/(kg ◦ C)−1

35

◦

C

Continued on next page 65

Page 65 of 72

Table A.1: Parameters and initial values of the Williams-Otto semi-batch reactor model.

value unit

1071

MA

molar mass A

100 kg/kmol

MB

molar mass B

100 kg/kmol

MC

molar mass C

ME

molar mass E

MG

molar mass G

MP

molar mass P

ρA,0

partial density A at t = 0

1 kg/m3

ρB,0

partial density B at t = 0

0 kg/m3

ρC,0

partial density C at t = 0

0 kg/m3

ρE,0

partial density E at t = 0

0 kg/m3

ρG,0

partial density G at t = 0

0 kg/m3

ρP,0

partial density P at t = 0

0 kg/m3

us

cr

heat transfer coeﬃcient

te

2.3082E-1 kJ/(m2 ◦ C s)−1

U

200 kg/kmol

an

200 kg/kmol 300 kg/kmol

d

M

100 kg/kmol

Ac ce p 1070

ip t

parameter description

reactor temperature at t = 0

V0

reactor volume at t = 0

Fscal

feed rate scaling factor

ρFA

density of feed rate FA (t)

1 kg/m3

ρFB

density of feed rate FB (t)

1 kg/m3

Appendix B

65

◦

TR,0

C

2 m3 1000 dm3 /m3

Algorithmic settings and initial values for numerical examples

66

Page 66 of 72

attribute value

integration tolerance 10−10

us

NLP solver SNOPT

cr

integrator S-LIMEX

ip t

Table B.1: Algorithmic parameters and tolerances for both illustrative example.

optimization tolerance 10−7

an

feasibility tolerance 10−6 max. ref. steps κmax , λmax

insertion tolerance ϵi

M

elimination tolerance ϵe

8

0.9

10−8

1

hor. refinement depth δh

0

te

d

ver. refinement depth δv

Ac ce p

Table B.2: Stopping threshold and constraint feasibility for multi-stage WO reactor.

attribute value

stopping threshold ϵΦ

10−5

constr. feasibility ϵI

0.01

Table B.3: Stopping threshold and constraint feasibility for multiple-shooting example.

attribute value stopping threshold ϵΦ

10−8

constr. feasibility ϵI

10−4

67

Page 67 of 72

Table B.4: Initial values for decision variables of the multi-stage problem.

k=1

k=2

ip t

phy. stage

1

num. intervals

16

1 1

us

num. stages l

cr

FA

5.784

0

approx. type

linear

constant

an

initial val. [dm3 /s ]

opt. variable

yes

no

num. stages l

1

1

num. intervals

1

16

initial val. [dm3 /s ]

0

5.784

approx. type

constant

linear

opt. variable

no

yes

num. stages l

1

1

num. intervals

16

16

initial val. [◦ C]

60.0

60.0

approx. type

linear

linear

opt. variable

yes

yes

dur. value 68 [s]

100.0

900.0

opt. variable

yes

yes

Ac ce p

te

d

M

FB

TW

∆tk

Page 68 of 72

ip t

k=1

FB 1

num. intervals

1

initial val. [dm3 /s]

1

1

1

1

5.784

5.784

5.784

yes

yes

yes

num. stages l

1

1

1

1

num. intervals

1

1

1

1

initial val. [◦ C ]

100.0

100.0

100.0

100.0

Ac ce p

te

TW

1

yes

d

opt. variable

k=4

linear linear linear linear

M

approx. type

5.784

k=3

1

an

num. stages l

k=2

us

phy. stage

cr

Table B.5: Initial values for decision variables of the multiple-shooting problem.

approx. type

linear linear linear linear

opt. variable

yes

yes

yes

yes

dur. value [s]

250.0

250.0

250.0

250.0

opt. variable

no

no

no

no

∆tk

69

Page 69 of 72

Evolution of identified control switching structure

ip t

for numerical examples Table C.1: Evolution of identified control switching structure for the multi-stage problem

cr

for the WO batch reactor. U: active upper bound, L: active lower bound, N: control

phy. stage

k=1

λ=1

FA

k=2

TW

U,L,N,U,L

TW

N,L,N,U N,L U,L,N U,L,U

N,L U,L,N

λ = 3 U,N,L,N,U,L

U,L,U

N,L U,L,N

d

M

λ = 2 U,N,L,N,U,L

FB

an

ref. step.

us

sensitivity-seeking or path-constrained arc.

κ = 1 U,N,L,N,U,L

U,L,U

N,L U,L,N

te

1073

Appendix C

κ = 2 U,N,L,N,U,L

U,L,U

N,L U,L,N

κ = 3 U,N,L,N,U,L

U,L,U

N,L U,L,N

κ = 4 U,N,L,N,U,L

U,L,U

N,L U,L,N

Ac ce p

1072

70

Page 70 of 72

ip t cr

Table C.2: Evolution of identified control switching structure for the multiple-shooting

control sensitivity-seeking or path-constrained arc.

ref. step. FB

k=2

TW

FB

k=3

an

k=1

TW

k=4

FB

TW

FB

TW

N

N,L

L

N,U

N,L,U N,L N,L

L

L,N,U

M

phy .stage

us

problem for the WO batch reactor. U: active upper bound, L: active lower bound, N:

U

N

U

L,U

λ=2

U

N,U,N

U

λ=3

U

L,U,N

U

N,U

U,L U,L

L

L,U

λ=4

U

L,U,N

U

N,U

U,L U,L

L

L,U

Ac ce p

te

d

λ=1

κ=1

U

L,U,N

U

N,U

U,L U,L

L

L,U

κ=2

U

L,U,N

U

N,U

U,L U,L

L

L,U

κ=3

U

L,U,N

U

N,U

U,L U,L

L

L,U

κ=4

U

L,U,N

U

N,U

U,L U,L

L

L,U

κ=5

U

L,U,N

U

N,U

U,L U,L

L

L,U

71

Page 71 of 72

*Highlights

ce pt

ed

M

an

us

cr

ip t

Exploitation of the control switching structure for multi-stage problems Exploitation of the control switching structure for multiple-shooting problems Wavelet-based grid refinement of nested multi-stage / multiple-shooting problem A new measure to evaluate the intermediate constraint violation

Ac

   

Page 72 of 72

Exploitation of the control switching structure in multi-stage optimal control problems by adaptive shooting methods

Exploitation of the control switching structure in multi-stage optimal control problems by adaptive shooting methods

Recommend Documents