Computational Algorithm of AOCP

For solving the optimal control problem, the essential idea of AOCP is to treat optimal control problems as initial-value problems by using iterative methods of nonlinear programming. The SQP method is selected to solve the nonlinear programming problems

transcribed from the discretization model of the original optimal control problem. Because SQP is a generalized gradient-descent optimization method and subsequently converges to a local rather than a global optimum, it solves the subproblem by providing both the direction of design improvement and the step size along the search direction. In this dissertation, the algorithms of AOCP and SQP are combined to form a general purpose solver, the OCP solver.

The architectural framework of the OCP solver, as illustrated in Figure 3.5, is composed of two computational blocks: the SQP algorithm and the OCP solver. Because the SQP algorithm is a well-known algorithm for optimization (Arora, 1989; Chong and Zak, 1996;

Rao, 1996), its implementation details can be found in a wide body of research and are therefore skipped in this dissertation. Basically, in each iteration of SQP, the values of the design variables are handed over to the OCP solver, which then uses them to calculate the values of the cost functional and the constraints. As shown in Figure 3.5, the OCP solver contains three major computational modules: discretization, calculation of current values of the state variables by applying the ODE solver, and estimation of the values of the cost functional and the constraints. Hence, the AOCP algorithm based on the SQP method can be described as following:

AOCP Algorithm:

Step 1. Choose b0, u0, Ns (maximum number of iteration)

N (number of time intervals)

ε

i (for convergence and stopping), k = 0 (iteration counter).

Step 2. Execute the discretization and parameterization process and calculate values for the following variables:

Time intervals T = [ T1, T2,…, TN ]^T defined in Eq. (3.10).

Interpolation parameters S^(k) by applying

u T

_k( )=

I S

( ^{( )k} , , )

T t

.

Design variable vector

P

^{( )k} = ⎣⎡

b

₁^{( )k} , ,

b

_π^{( )k} ,

T

₁^{( )k} ,

T

_N^{( )k} ,

S

₁^{( )k} , ,

S

_m^{( )k}_N⎤⎦ . ^T

Step 3. Determinate the state variables x^(k) by solving the initial-value problem defined in Eqs. (3.16) and (3.17) with P^P(k).

Step 4. Calculate the values of the cost functional Eq. (3.13) and the constraints Eqs.(3.14) and (3.15).

Step 5. Calculate the gradients of the cost functional and the constraints.

Step 6. Find the descent direction d by solving the QP subproblem defined in Eqs. (3.6) and (3.7).

Step 7. Check feasible and convergence criteria.

(a) Convergence for SQP:

IF (KT condition is satisfied) THEN

Algorithm converged, Stop.

(b) Stopping criteria:

ΔP = P^P(k+1) –P^(k)P

IF (ΔP

^PTΔP ≤

ε

i ) THEN Stop. (design variable not changing)

IF (k = N

s) THEN Stop. (maximum iteration reached)

Continue

Step 8. Calculate the step size α^(k).

Step 9. Update the Hessian matrix H ^(k) by applying a BFGS approximation defined in Eq.

(3.8).

Step 10. k = k +1; Go to Step 2.

3.5 Summary

This chapter has introduced the numerical preliminaries – including NLP formulation, the SQP method, the control parameterization technique, and dynamic constraint treatments – for developing the OCP solver. Also discussed were methods for solving the NLP and the

first-order necessary condition for continuous constrained NLP problems. The text also introduced the SQP method that serves as the kernel of proposed method, and provided its algorithm. In addition, because this dissertation uses control parameterization with an interpolation function to decrease the numbers of design variables and so make the solver more efficient, discretization and parameterization techniques that transcribe the optimal control problems into NLP problems were developed. Also introduced were several numerical schemes involved in the proposed solver – including ODE solvers and integration and interpolation schemes – and finally, the computational algorithm of the ACOP method that will help implement the OCP solver.

Figure 3.1 Methods for continuous constrained NLPs (Wu, 2000).

Figure 3.2 Conceptual flowchart of the SQP method.

Initialization

Initial guess x0

Hessian matrix H = I Iterative index k = 1

Calculate values of cost function and constraints.

Calculate gradients of cost function and constraints.

Calculate the descent direction, d^(k)

Meet the convergence criteria ?

Line search

Search step size, α^(k) Update Hessian matrix

H

^(k)

Update design variables x^(k)

x

^(k+1) = x^(k) + α^(k)

d

^(k)

k = k + 1

Show results Yes No

Figure 3.3 Problem-transcribing Process.

Φ

t

Conventional formulation

Φ

t

Worst-case design formulation

Φ

t

Hybrid formulation

Figure 3.4 Dynamic constraint treatments.

Initialization

P⁽⁰⁾, H⁽⁰⁾= I, k = 0

Gradients information

Descent direction convergence?

Searching step size Update Hessian matrix Update design variables

Show results No

k = k + 1

Solve IVP

Performance index J0* Constraints J_i*

Sequential Discretization

X^(k) P^(k)

Yes

(SQP) AOCP

Figure 3.5 Conceptual flowchart of the AOCP method.

CHAPTER 4

A CONVENIENT SOLVER FOR SOLVING OPTIMAL CONTROL PROBLEMS

4.1 Introduction

Even though, over the last two decades, theoretical and numerical methods for solving optimal control problems have been extensively studied and many well-designed algorithms have been proposed, engineers must still expend much effort to reformulate the nonlinear programming problems for different control problems. On the other hand, reworking the corresponding programs for the nonlinear programming problem is a waste of time and even more tedious. Therefore, developing a general OCP solver that offers a systematic process for solving various optimal control problems has become imperative for engineers, particularly for those who are inexperienced in optimal control theory or numerical techniques.

As mentioned in the previous chapter, many well-developed subroutines or program units for numerical analysis are involved in developing a general OCP solver, e.g., integration routines, ODE solvers, interpolation schemes, and general constrained optimization solvers. A general OCP solver, named the OCP solver, was designed that consists of the subsystems or components for such numerical subroutines. Its modular programming features enable the OCP solver to choose different numerical schemes flexibly and be easily upgraded by replacing subroutines with new versions. The implementation details of the kernel module and user interface of the OCP solver are presented in following sections.

4.2 Multifunctional Optimization System Tool - MOST

In the AOCP method, the optimal control problem is converted into an NLP program so that any reliable nonlinear constrained optimization solver can be applied to solve it numerically. A great deal of attention has been paid to using the SQP method to solve NLP problems (Tseng, 1987; Jaddu and Shimemura, 1999). For this dissertation, the

multifunctional optimization system tool MOST (Tseng, 1993), based on the SQP method, has been chosen to solve the NLP problem. This powerful optimization software was developed to solve multi-objective optimization problems with both continuous and discrete design variables (Tseng et. al., 1993). The MOST software contains three main modules for dealing with continuous variables, discrete variables, and multi-objective optimization, respectively.

In the primary module, a SQP method (Arora, 2004) is employed to perform the single-objective optimization for problems with continuous design variables. The SQP is selected because of its accuracy, efficiency, and robustness. MOST’s accuracy and stability has been tested in the research using 115 test problems with 2 to 96 design variables given by Hock and Schittkowski (1980) The results were satisfactory and also indicated that MOST can handle large-scale engineering optimization with excellent convergence (Tseng et. al., 1988a; 1988b). To cope with the discrete-valued optimal control problems arising from discrete design variables, an enhanced branch-and-bound method (Tseng et. al., 1995) was integrated into the program. In this module, the original design space of discrete variables is converted into one with continuous variables by dropping the noncontinuous restrictions sequentially. In each of the converted continuous design spaces, the SQP module described above is then utilized to find the optimal values.

Nevertheless, in many engineering applications there frequently exist several mutually conflicting or competing objectives and requirements. Therefore, multi-objective (vector) optimization offers a very promising way to handle such problems. For multi-objective optimization, MOST provides decision makers, goal programming, compromise programming, and the surrogate worth trade-off method (Evans, 1984; Tseng and Lu, 1990) to help users determine the best compromised solutions to nonlinear problems.

It is known that a rigorous formulation of the design problem helps the designer better understand the problem and a proper mathematical formulation leads to a good solution.

MOST provides an input data file that includes the initial design, and it transcribes the design problem by coding subroutines that include evaluation of the cost, in routine cusermf, and constraint functions, in routine cusercf (see Figure 4.1). Because of its user-defined subroutines, MOST can be extended to flexibly integrate other analysis packages, e.g., ANSYS (Yang et al. 1992; Lin et al., 1992), EUCLIS-IS (Wang, 1993), and MATHEMATICA (Su, 1994). Moreover, with the assistance of an interface coupler, MOST can deal with the complexity and large size of engineering systems that have no explicit relationships between inputs and system outputs (Huang, 1994). The architecture of the MOST interface coupler is shown in Figure 4.2. The IAOS, a new distributed version of the interface coupler, was developed to deal with analysis packages installed on different machines (Huang, 1994). In IAOS, MOST and the interface coupler are merged into a powerful new optimizer (see Figure 4.3).

4.3 Structure of the Proposed OCP Solver

The kernel of the OCP solver is written in FORTRAN and has been tested on a UNIX platform. Because of the increased popularity and technical developments in the calculation ability of the personal computer (PC), for this dissertation, the OCP solver has been transplanted onto a PC platform. The structured chart for the entire OCP solver is given in Figure 4.4, and the connections between the OCP solver and MOST are shown in Figure 4.9.

Four primary independent modules of the OCP solver are introduced below.

CTRLMF module: The structure of this module, given in Figure 4.5, links with MOST’s

user-defined subroutine cusermf and contains six subroutines. The control routine CTRLMF is used to calculate the value of the performance index, to partition the array, and to check available memory. The pseudocode for the CTRLMF module of the AOCP algorithm is shown in Table 4.1.

CTRLCF module: This module, whose structure is shown in Figure 4.6, calculates the

constraint function values. First, the information on the state trajectory is passed from the

CTRLMF module and then the values of the functional constraints and dynamic constraints

are calculated in routine CTRLCF. In this module, subroutine PTCST is used to calculate the number and values of dynamic constraints for the alternative treatments described in Section 3.4.2.

CTRLMG module: This module, illustrated in Figure 4.7, calculates the design derivatives of

the performance index by using the direct differentiation or adjoint variable method. In the

DDMMG routine, the design derivatives are calculated using direct differentiation (DDM); in

the AVMMG routine, by the adjoint variable method (AVM). The ANSAVM routine determines the terminal conditions and integrates the adjoint differential equation using a backward numerical integration scheme.

CTRLCG module: This module, outlined in Figure 4.8, is in charge of design sensitivity

analysis, which determines the effect of a change in the current design on the performance index and constraint functions. Two common methods, DDM and AVM, described in Section 3.4.3, are implemented in this module. The control routine DDMCG is evaluated to obtain the design derivatives using DDM. The information for state variable derivatives with respect to the design variables is then passed from the CTRLMG module. If the AVM is needed to calculate the design derivatives, the AVMCG routine is executed.

In addition to these four modules, many useful routines are shared by a variety of modules, e.g., DIFSOL, INGSOL, DDERKF, DDEABM, DDEBDF, SIMPSN, GAUSS, TBFIT,

TGVAL, and CTRLF. The routine DIFSOL contains three differential equation solvers (RKF,

ABM and DBF), and the routine INGSOL is used to calculate the value of an integral. In

INGSOL, the SIMPSN (bases on Simpson’s rule) and GAUSS (following the Gaussian

quadrature formula) are called on to evaluate the integral. Both the DDERKF and the

DDEABM can integrate a system of first-order differential equations, the first using the

Runge-Kutta-Fehlberg method and the second using the Adams-Bashforth-Moulton predicator-corrector formulas for orders one through. The latter implements a backward differentiation formula in the routine. The TBFIT is used to calculate the coefficients of the interpolation function and the TGVAL, to calculate the value of functions and their derivatives.

4.4 The OCP Solver in Cooperation with MOST

Because of the extension flexibility of MOST, the OCP solver is herein treated as a MOST module. The linkage of MOST to the OCP solver is composed of four user-defined MOST subroutines: cusermf, cusermg, cusercf and cusercg. The OCP solver also has four user-defined subroutines for defining an optimal control problem: USERMF, USERMG,

USERCF and USERCG. The architecture for MOST and the OCP solver is illustrated in

Figure 4.9.

In the OCP solver, the subroutine USERMF provides the cost functional value (performance index value); subroutine USERCF provides the constraint function values; and subroutines USERMG and USERCG provide the cost function gradient and the gradients of the active constraints, respectively. A fifth subroutine, USEROU, can be developed by the user to perform a subsequent optimality analysis for the optimal solution and obtain more output.

If the analytic expressions for the gradients in USERMG and USERCG are not available, MOST provides the option of calculating gradients using a finite difference method that can be specified as forward, backward, or central. More details are provided in the MOST 1.1 User’s Manual (Tseng et al., 1993). As the architecture given in Figure 4.9 shows, connections exist between the optimizer MOST, the five user subroutines, and the four independent modules, CTRLMF, CTRLMG, CTRLCF and CTRLCG, for the OCP solver. The four modules contain the performance index, the functional and dynamic constraints, the gradient of the performance index, and the constraint function gradients for the optimal control problems, respectively. They can be connected to each other to form a general

constrained optimization solver.

4.5 User Interface for the OCP Solver

The user interface for the OCP solver consists of two parameter files and four subroutines.

Users can specify the optimization parameters and numerical schemes in the parameter files of both MOST and the OCP solver, such as acceptable violation of constraints for feasible designs, the differential equation solver, the integration rules, and the interpolation scheme.

The user interfaces for the OCP solver and MOST are shown in Figure 4.10.

The MOST optimal parameter file, shown in Table 4.2, is used to configure the optimizer parameters (details can again be found in the MOST user manual [Tseng et. al., 1993]).

Details of the OCP parameter file fort.11 – which contains information on the number of equations, grid points, equality and inequality functional constraints, and the parameters for numerical schemes – can be found in Tseng (1987). Table 3 gives the OCP parameter file for the van der Pol oscillator problem to be introduced in Section 4.7.1. Once all four necessary user-defined subroutines, FFN, GFN, HFN, and Z0FN, are ready, they should be linked to the OCP through the MOST kernel. Figure 4.10 shows the relationship between the user-defined routines and the MOST modules. FFN evaluates the integral terms of the performance index or functional constraints, while GFN calculates the values of the first term of the performance index or functional constraints and the dynamic constraints. HFN is used to evaluates the system state equation f(b, u(t), x(t),t), and finally, Z0FN is responsible for calculating the initial values of the state variables, x(t0). Figure 4.11 gives a flowchart for the OCP solver.

4.6 Systematic Procedure for Solving the OCP

In this dissertation, the OCP is converted into an NLP problem using an admissible optimal control problem formulation then the optimizer based on the SQP method is used to solve the NLP problem numerically. These procedures can now be directly implemented, and

the complicated details of transformation and programming automatically completed, in the proposed OCP solver. Because the optimal control and state trajectories are obtained and recorded in the output files, engineers can follow an efficient and systematic procedure to solve various optimal control problems. The procedure for solving the OCP with the OCP solver is as follows:

1) Defining the OCP problem following the formulation defined in Section 2.2.

2) Preparing the parameter files and user-defined subroutines according to the formulation.

3) Compiling the user’s subroutines and linking with the OCP solver.

4) Executing the OCP solver and obtaining the optimal results.

4.7 Illustrative Examples

Two types of optimal control problems mentioned in the literature have been used as test problems to evaluate the performance of the proposed method. In the AOCP method, both the acceptable violation of constraints for feasible designs and the acceptable tolerance for the convergence parameter are 10^-3. The numerical results for all example problems were obtained on a Pentium 4 Celeron 1.2 GHz computer with 384 MB of RAM.

4.7.1 The van der Pol Oscillator Problem

The van der Pol oscillator problem was given and solved by Bullock and Franklin (1967) using a second variation method. The problem was also used by Jaddu and Shimemura (1999) to verify their computational method. In this dissertation, it is further used to evaluate the performance and capabilities of the proposed method and the OCP solver. The van der Pol oscillator problem can be formulated by the following minimization

dt

subject to

Based on this problem, Jaddu and Shimemura considered three cases that can also be solved by the OCP solver: the unconstrained problem, the terminal state constrained problem, and the terminal states and control constrained problem.

Case I: Free end point and no control constraints

The optimal solution for this problem found by Bullock and Franklin (1967) using a second variation method was J0* = 1.433508, while that found by Jaddu and Shimemura (1999) using a ninth-order Chebyshev series to approximate x1(t) was J0* = 1.4334872.

Using the OCP solver, in which the control variable u is discretized into 21 grid points, the optimal value is J0* = 1.4334723, smaller than both earlier reported results. The numerical parameters for MOST are listed in Table 4.2, and the optimal control and state trajectories are shown in Figure 4.12.

Case II: Terminal state constraint

( ^x ( ) ^t

^f

) ^{1 x}

²

( ) ^t

^f

^x

¹

( ) ^t

^f

ϕ

= − + = 0

(4.3)

For this problem, Bullock and Franklin (1967), again using the second variation method, found an optimal value of J0* = 1.6905756, while Jaddu and Shimemura (1999), also using a ninth-order Chebyshev series to approximate x1(t), found an optimal value of J0* = 1.6857113. In this study, the terminal state constraint is treated as an equality constraint and the other number parameters, the same as in case I. With the OCP package, the value obtained is J0* = 1.6856957. Figure 4.13 shows the optimal control and state trajectories for the proposed OCP solver.

Case III: Terminal state constraints and saturation constraints on control

The terminal state constraints and the saturation constraints on control are described in the following equation:

When Bashein and Enns (1972) solved the problem, they obtained J0* = 2.1439039, while Jaddu and Shimemura (1999), this time using a twelfth-order Chebyshev series to approximate x1(τ), found an optimal value of J0* = 2.1443893. The solution produced by the OCP solver is an optimal value of J0* = 2.1375360. The optimal control and state trajectories for the OCP solver are shown in Figure 4.14.

4.7.2 Time-optimal Control Problem: Overhead Crane System

Overhead cranes are widely used in factories and workplaces to transport objects. An overhead crane system, like that sketched in Figure 4.15, is a high-order nonlinear system that consists of a cart with a point load suspended by cables. The control problem is to transfer the load from an arbitrary point A to point B in minimal time subject to the requirement of zero residual vibration at point B. The control inputs are the horizontal acceleration of the cart and the hoisting acceleration of the cable. Hu et al. (2002) proposed this problem and solved it using an enhanced DCNLP method. In this dissertation, this problem will be used to demonstrate the ability of the proposed method to solve a high-order time-optimal control problem.

Given x₁ =z, x₂ =z,

x

₃ =

θ

,

x

₄ =

θ

=

ω

,

x

₅ =

l

,

x

₆ =

l

as the state variables, and u₁ = z, as the control inputs, the OCP formulation of the overhead crane system can be minimized as follows

l

u

₂ =

J0 = tf (4.5) 4, 0]^T where g is the gravitational acceleration.

The state and control constraints are as follows:

tf

Using the admissible control formulation delineated in the previous chapter, the control variables are converted into design variables that can then be treated as design variable boundaries. Furthermore, the state constraints are transferred into standard constraint form as follows:

In this problem, both the state and the control variables are divided into 101 grid points.

The minimum time J* = tf = 12.0004 is solved in the OCP solver by applying a cubic piecewise interpolation scheme to the control function. Two local optimal solutions are obtained by the OCP solver with different initial points. The trajectories of the rope angle and angular velocity are shown in Figure 4.16, in which the solid line represents the results obtained by the DCNLP method (Hu. et al., 2002) and the dashed line represents a second

optimal solution obtained by the OCP solver. As the figure illustrates, the solid line totally matches the results obtained by Hu et al. (2002), meaning that one of the optimal solutions found by the OCP solver tallies exactly with the trajectories obtained by the DCNLP method (Hu. et al., 2002). In addition, the performance index (terminal time, tf) obtained by the OCP solver (tf ≅ 12.00) is very close to the result using the DCNLP method (tf ≅ 12.00). Moreover, according to the trajectories shown in Figure 4.16, the amplitudes of rope angle and angular velocity for the second optimal solution obtained by the OCP solver, as represented by the dashed line, are smaller than the others. Figures 4.17 and Figure 4.18 depict the corresponding inputs and states with local optimal solutions, respectively. In Figure 4.17, the trajectories of the control inputs conform to the dynamic control constraints given in Eq. (4.8). According to

optimal solution obtained by the OCP solver. As the figure illustrates, the solid line totally matches the results obtained by Hu et al. (2002), meaning that one of the optimal solutions found by the OCP solver tallies exactly with the trajectories obtained by the DCNLP method (Hu. et al., 2002). In addition, the performance index (terminal time, tf) obtained by the OCP solver (tf ≅ 12.00) is very close to the result using the DCNLP method (tf ≅ 12.00). Moreover, according to the trajectories shown in Figure 4.16, the amplitudes of rope angle and angular velocity for the second optimal solution obtained by the OCP solver, as represented by the dashed line, are smaller than the others. Figures 4.17 and Figure 4.18 depict the corresponding inputs and states with local optimal solutions, respectively. In Figure 4.17, the trajectories of the control inputs conform to the dynamic control constraints given in Eq. (4.8). According to

在文檔中 動態系統最佳化設計之計算方案與其應用 (頁 54-0)