• Author / Uploaded
• Utku Oğuzman

Citation preview

ENGINEERING TOOLS, TECHNIQUES AND TABLES

OPTIMAL CONTROL ENGINEERING WITH MATLAB

No part of this digital document may be reproduced, stored in a retrieval system or transmitted in any form or by any means. The publisher has taken reasonable care in the preparation of this digital document, but makes no expressed or implied warranty of any kind and assumes no responsibility for any errors or omissions. No liability is assumed for incidental or consequential damages in connection with or arising out of information contained herein. This digital document is sold with the clear understanding that the publisher is not engaged in rendering legal, medical or any other professional services.

ENGINEERING TOOLS, TECHNIQUES AND TABLES Additional books in this series can be found on Nova‟s website under the Series tab.

Additional e-books in this series can be found on Nova‟s website under the e-book tab.

ENGINEERING TOOLS, TECHNIQUES AND TABLES

OPTIMAL CONTROL ENGINEERING WITH MATLAB

RAMI A. MAHER

New York

Copyright © 2013 by Nova Science Publishers, Inc. All rights reserved. No part of this book may be reproduced, stored in a retrieval system or transmitted in any form or by any means: electronic, electrostatic, magnetic, tape, mechanical photocopying, recording or otherwise without the written permission of the Publisher. For permission to use material from this book please contact us: Telephone 631-231-7269; Fax 631-231-8175 Web Site: http://www.novapublishers.com

NOTICE TO THE READER The Publisher has taken reasonable care in the preparation of this book, but makes no expressed or implied warranty of any kind and assumes no responsibility for any errors or omissions. No liability is assumed for incidental or consequential damages in connection with or arising out of information contained in this book. The Publisher shall not be liable for any special, consequential, or exemplary damages resulting, in whole or in part, from the readers‟ use of, or reliance upon, this material. Any parts of this book based on government reports are so indicated and copyright is claimed for those parts to the extent applicable to compilations of such works. Independent verification should be sought for any data, advice or recommendations contained in this book. In addition, no responsibility is assumed by the publisher for any injury and/or damage to persons or property arising from any methods, products, instructions, ideas or otherwise contained in this publication. This publication is designed to provide accurate and authoritative information with regard to the subject matter covered herein. It is sold with the clear understanding that the Publisher is not engaged in rendering legal or any other professional services. If legal or any other expert assistance is required, the services of a competent person should be sought. FROM A DECLARATION OF PARTICIPANTS JOINTLY ADOPTED BY A COMMITTEE OF THE AMERICAN BAR ASSOCIATION AND A COMMITTEE OF PUBLISHERS. Additional color graphics may be available in the e-book version of this book.

Library of Congress Cataloging-in-Publication Data ISBN:  (eBook)

Published by Nova Science Publishers, Inc. † New York

To my students over the years

CONTENTS Preface

xi

Acknowledgments

xiii

Chapter 1

Mathematical Background and Optimal Problem Modeling 1.1. Introduction 1.2. Optimal Control in Engineering 1.3. Controller Design Concepts 1.4. Finite Dimensional Optimization–A Mathematical Review 1.5. Optimal Problem Modeling Problems References

1 1 2 3 4 9 30 32

Chapter 2

Controller Design Based on Parameter Optimization 2.1. Introduction 2.2. Performance Indices 2.3. Preliminary Design Concept – The Inward Approach 2.4. Parameter Optimization Design – Outward Approach 2.5. Limitations of Parameter Optimization 2.6. Parameter Optimization with Constraints 2.7. Control Vector Parameterization 2.8. Parameters Optimization via Genetic Algorithms Problems List of MATLAB Programs References

35 35 35 41 44 55 64 70 75 79 84 94

Chapter 3

Calculus of Variations 3.1. Introduction 3.2. Motivation 3.3. Standard Problem Statement 3.4. Transverzality Conditions 3.5. Extreme Functional with Dynamic Constraints Problems References

95 95 95 97 107 113 118 124

viii

Contents

Chapter 4

Optimal Control Based on Calculus of Variations 4.1. Introduction 4.2. Optimal Criteria 4.3. Statement of Optimal Control Problem 4.4. Necessary Conditions of Optimality 4.5. Terminal Controllers and Regulators 4.6. Terminal Constraints and Variable Terminal Time 4.7. Approximate Optimal State-Feedback Control 4.8. Numerical Solution Using MATLAB TBVP Solver Problems List of MATLAB Programs References

125 125 125 127 132 145 155 164 167 176 183 189

Chapter 5

Optimal Control with Input and State Variable Constraints 5.1. Introduction 5.2. Input Constraints Optimization 5.3. Pontryagin’s Principle 5.4. Optimal Time Problem 5.5. Minimum Control-Effort Problems 5.6. Singular Problem 5.7. Inequality Constraint of Pure State Variable Problems List of MATLAB Programs References

191 191 191 193 197 211 230 242 250 255 258

Chapter 6

Dynamic Programming 6.1. Introduction 6.2. Multi-Stage Decision Problem and Principle of Optimality 6.3. A Simple Control Problem in Discrete-Time 6.4. Continuous Form of Dynamic Programming 6.5. Parametric Expansion 6.6. Cross-Product Terms 6.7. Differential Dynamic Programming DDP Problems List of MATLAB Programs References

259 259 259 261 269 280 289 290 300 306 313

Chapter 7

Linear-Quadratic (LQ) Optimal Control 7.1. Introduction 7.2. Types of Optimal Control Problems 7.3. Optimal Solution of the State Regulating Problem 7.4. Selection of the Q and R Matrices 7.5. Optimal Solution of the Output Regulating Problem 7.6. Optimal Solution of the Tracking Problem 7.7. Equivalent Quadratic Cost Functional 7.8. Solution Based on Hamiltonian Matrix for Invariant Systems 7.9 Discrete-Time LQ Optimal Control

315 315 315 316 323 328 335 348 353 360

Contents

Chapter 8

Index

ix

7.10. Iterative LQ for Nonlinear Systems Problems List of MATLAB Programs References

367 372 379 385

Optimal Solution Based on Genetic Programming 8.1. Introduction 8.2. Genetic Programming Review 8.3. Enhanced Genetic Programming 8.4. Synthesizing Input Constraints Optimal Control Problem Problems References

387 387 387 391 404 416 417 419

PREFACE People aspire after the optimum when working, or doing tasks, or completing duties. However, writing a book about optimal control for engineering students would ask for careful treatment. This is deliberately because the subject is so wide-ranging and technically enriched and new branches are still emerging in this field. Furthermore, the optimal theory always finds new areas of challenges in different fields. For instance, the motor control is one of the newest areas in the bio-field. As I lectured the optimal control subject over the years, graduate engineering students often asked the same questions: What is the difference between an optimized controller and an optimal controller? How can engineers implement optimal controllers? Is there a bridge between conventional control design methods and the optimal control approaches? Is it worth learning this subject first, before studying the robust optimal control approaches? How does studying this subject help with the final project or with future control engineering practice? How the benefit of the MATLAB facilities for solving optimal controller design problems can be exploited? These questions made me think of writing a book addressed to graduate engineering students to satisfy their engineering inquisitiveness. This book has been written primarily for use as a first-year-graduate-level optimal control course. It covers the deterministic part of the optimal control theory. Essential theories and explanations of the various techniques are explained in a simple manner without rigorous mathematics but also without sacrificing accuracy. A sufficient number of worked examples and problems are included to highlight the engineering design concepts. Many of these examples are worked by utilizing simple and readable MATLAB programs. Simulink models are also used on one hand to provide and check solutions, and on the other hand, to ensure understanding of the implementation aspects. The book consists of eight chapters. Chapter 1 introduces optimal control principles, essential definitions, and concepts of optimal controller design, and a mathematical review of static optimization. The modeling of optimal control problems for some engineering applications is finally presented. In chapter 2, the design of controllers based on parameter optimization via the inward and outward approaches is considered. The performance index as a measure of the quality of the required response and the resultant forms of the closed-loop are first considered. While the inward approach is briefly introduced, the outward is thoroughly discussed, and many illustrative examples are given. Control vector parameterization and genetic algorithm are also presented. To introduce the optimal control theory from a mathematical point of view, Chapter 3 reviews the approach of calculus of

xii

Rami A. Maher

variations. The reader will sense the close relation between the material of the chapter and the optimal control framework presented in the next chapters. Chapter 4 is intended to highlight the design concepts of optimal control based on calculus of variations. The necessary and sufficient conditions of optimality of various types of problems are given and several examples, which are illustrated with the aid of MATLAB programming, are presented. Furthermore, the chapter includes the approximate technique of optimal or suboptimal state control law and the numerical solution of the Two-Point Boundary Value Problem (TBVP) using the available MATLAB special functions. The very famous optimality principle of Pontryagin is presented in chapter 5, which covers the problem of constrained input or state imposed on the control systems. It includes the solutions of optimal time problems, the minimum control-effort problems, singular problems, and the problem of state variable inequality constraints. Analytic optimal solutions are presented, as well as the numerical gradient method that is applied to solve the state variable inequality problem. Chapter 6 introduces the dynamic programming approach for both of discrete and continuous systems, where the “principle of optimality” introduced by Bellman is thoroughly discussed and illustrated. The parametric expansion approach is considered as a possible tool for solving the Hamilton-Jacobi-Bellman (HJB) equation (the necessary condition of optimality). In addition, differential dynamic programming (DDP) is included as an approach for nonlinear systems. Nevertheless, the results for the Linear Quadratic Regulator (LQR) approach are derived via the application of the optimality principle and the topic of linear quadratic optimal control is found in chapter 7. In chapter 7, the mathematical formulation of three basic control problems is addressed: the regulating problem, the output-regulating problem, and the tracking problem. The stability of the optimal solution, equivalent cost functions, solutions based on the Hamiltonian matrix are all considered in this chapter for both continuous-time and discrete-time LQR problems. The chapter closes with a new iterative linear quadratic regulator (ILQR) approach adopted for nonlinear systems. Chapter 8 covers the optimal control solution based on one of the newest evolutionary techniques, the genetic programming (GP). The chapter includes a brief introduction for GP and details of an enhanced algorithm to solve the optimal problem. Both input unconstrained and input constrained optimal problems are considered. Rami A. Maher

ACKNOWLEDGMENTS Thanks for God I wish to thank my merciful wife, Basma, great son, Ali, and lovely daughters, Raya and Leena, for their patience, understanding, and encouragement in the long preparing of the manuscript. Special dedication goes to Professor Dr. Mike Johnson, who gave me worthy recommendations and supports. Also I would like to thank dear friend Professor Dr. Waladin Said for his help during the early preparing of the book manuscript.

Chapter 1

MATHEMATICAL BACKGROUND AND OPTIMAL PROBLEM MODELING 1.1. INTRODUCTION One of the engineering challenges is to perform tasks in a best possible way with a small amount of energy and with a minimum number of resources. Nature of things makes the meaning of the word "best" very ambiguous and confusing. Moreover, tasks spread widely from the manufacturing of simple objects to complex systems. For instance, to obtain a maximum number of specified circle pieces out of a given rectangular sheet is a simple example. On the other hand, to control complicated dynamic systems like minimizing the miss distance between a pursuer missile and its target is a complex one. Therefore, the way, or approach to undertake one task is not necessarily preferable for other tasks. However, for a specific system or process, a precise definition of the demands, the available resources, and the physical constraints make engineers more comfortable and hopeful about finding the best, most favorable, most advantageous or optimal solution. The practical realization of this optimal solution is the next stage the engineer has to accomplish. During implementation, engineer encounters often many obstacles that oblige him to postulate a near or sub-optimal solution. It can be reached by modifying or discarding some of the customer demands without much loss. The extreme (minimum or maximum) of a goal function, a cost function, performance index, or satisfying certain optimal criterion is the optimization process. Intuitively, a maximization problem can be put easily as minimization one and vice versa. The alternative usage of terminology to express almost the same meaning is due to the wide range of applications. In this book, the optimization process considers the problem of finding the optimal input control for a dynamic system provided that a specific optimal criterion and a set of constraints are satisfied. Historically, the calculus of variations has constituted the theoretical framework for solving the optimal control problems. The most important necessary and sufficient conditions were derived based on the concepts and theories of this branch. Later, several important conductive results and new approaches to tackle the problem were established. The Pontryagin‟s principle of maximum and Bellman principle of optimality are the most famous.

2

Optimal Control Engineering with MATLAB

Nowadays, the optimal robust control concepts and theories make a burst through in the field where multiple objectives can be handled simultaneously. To pass over these entire subjects, it is believed that for better understanding one should start from the beginning. For this reason in this chapter some fundamental concepts, essentials of optimization theory and optimal problem modeling are presented. Chapter two considers controllers design based on the methodology of parameter‟s optimization. Chapter three is devoted to explain briefly the fundamental concepts of the calculus of variations and how it relates to the optimal control problems. Chapter four explains and discusses solutions based on the calculus of variation methods. Chapters five and six respectively treats the problem based on the vision of Pontryagin and Bellman. The LQR technique, which is one of the most important results in the field, is considered in chapter seven. The last chapter covers partially one of the most recent evolutionary techniques, genetic programming GP. The optimal solution is proposed via an enhanced algorithm of GP.

1.2. OPTIMAL CONTROL IN ENGINEERING Optimal Control theory is a mature mathematical discipline with numerous applications. In fact, the applications of optimal control theory are very widespread that includes different fields. It is not exclusive to the engineering area. To give some idea of the breadth of this field, we shall make a few sample citations. For instance, the treatment of a two-strain tuberculosis model [1] or drug treatment strategies of the HIV immunology model [2] are two examples of many applications in the biological area. Besides engineering problems, mathematicians are usually interested in financial models, social models, production and maintenance models [3], transportation or vehicle-passenger models [4], etc. Mathematicians also contribute in many different theoretical new approaches and develop modern solution methods. For example, the solution of the famous Hamilton-Jacobin-Bellman HJB equation, the solution of optimal singular and bang-bang problems, the sensitivity of the optimal solution to the variation of model parameters, etc., has taken the greatest attention [5, 6, 7, 8]. Engineers believe that tuning the controller parameters by a trail-and-error procedure (a usual engineering practice) in some applications is a good motivation to search for an optimal solution. On the other hand, engineers do not look at the mathematical solution as an end of the story but what is more important to them is the possibility of implementing that solution in real time with the available resources. Therefore, engineers contribute positively in the developments of the optimal control field. Weighting of a control cost versus quality in the performance function is one of the fundamentals in optimal control theory, however, stating a new cost function that considers the cost of the sensors and actuators to design a special optimal controller [9] is an example of such engineering contributions. Obviously, it is hard to list or to reference all engineering applications, but it will be illustrative, to mention some abstract examples. Optimal temperature control of a drug store is an engineering task that belongs to biological fields. Suppose that the outside temperature is higher than the required inside temperature while a demanded electrical power to drive the cooling compressor can be chosen to range from minimum to maximum values. The store size, the outside temperature,

Mathematical Background and Optimal Problem Modeling

3

and the cooling system specifications state a dynamic system. Thus, the behavior of the inside temperature will depend on this system. The customer may wish to keep always the variation of the inside temperature as a minimum as possible without regarding the cost of the consumed electrical power. Alternatively, he wishes to make a tradeoff between the consuming power and the inside temperature increment but without overshooting the required value for a specified period of time. Control engineer will look to the former wish, as it is easier to fulfill. In modern applications, one can realize that the optimal control problem often shares with many engineering branches. Fuel-optimal control of continuously variable transmissions CVT power trains [10] is an example, where the solution of the fuel problem for transient conditions is assumed a challenging problem. Another engineering stimulating example is to bring a spacecraft to a soft landing on the lunar surface using the least amount of fuel, i.e. minimum fuel used up or also to maximize the amount remaining once it has landed. A comprehensive mathematical model with so many constraints describes such a problem. In the robotics field, one can easily find that the optimal control theories are extensively used to solve different problems of robot manipulator motion, real-time trajectory planning, force control, lateral and velocity control in car-like robot, etc. Most of these problems are defined by minimax optimal criteria, and they consider stochastically disturbances to design robust controllers [11]. The fact that the optimal control theory can be applied irrespective of the order and the nature of the plants pushes engineers to design and elaborate powerful robust controllers of a wide range of industrial applications [12]. For example, one real airplane model may consist of twelve first-order nonlinear differential equations and the task of designing a robust controller is a challenging task. Recent application of optimal control theory is extended to switching systems, which are a particular class of the hybrid system that consists of several subsystems and switching laws orchestrating the active subsystem at each instant of time. Many real-world processes such as chemical processes, automotive systems, and manufacturing processes, etc., can be modeled as switching systems. Optimal control problems for switched system, which require the solution of both optimal switching sequences and the optimal continuous inputs, have attracted many engineers and researchers [13].

1.3. CONTROLLER DESIGN CONCEPTS The fundamental problem of designing a controller is the attempt to fulfill a set of time and frequency specifications. They define or state the required performance of the control system in terms of certain measurable quantities. The rising or regulating time, settling time, dynamic and steady-state errors, maximum overshoot, phase and gain margins, resonant peak or maximum singular value, resonance frequency or system bandwidth, and the sensitivity characteristics, are the most usual specifications. It is a known fact that some of these quantities are contradicting in the sense that the response improvement cannot be carried out in all directions. In practice, it is always the case, that the customer requirements cannot be satisfied completely by a simple gain adjustment or even by a simple controller. Although most

4

Optimal Control Engineering with MATLAB

controlling schemes, which are based on classical formulas, provide a stable response, the designer may also suffer from the burden of trials to reach an accepted tradeoff between the contradicting time and frequency specifications. Moreover, due to the rapid success of technology in producing sensors and actuators, passive and active elements, microelectronics, the computation tools, etc., the demands‟ roof is as well raised. Therefore, the usual tradeoff approach, which satisfies some of the contradicting requirements, becomes less acceptable and at the same time, control engineers are encouraged to face challenging demands. In general, there are two approaches of designing a specific controller the classical openloop synthetic design, and the analytic closed-loop design. For the latter, the design is based on two approaches. The first approach is to select the configuration of the plant and the compensator interconnection, then to determine the parameters of the compensator, which meets the required specifications. The second one is to find that system which fulfils given specifications, and then the necessary compensator is computed. An alternative approach is to conform to the minimizing or maximizing a certain performance index without considering any specific configuration. This approach affirms the design of the optimal controllers. For optimal controller designs, there are also two approaches to handle the controller design problem. The first approach is based on adopting a fixed controller structure as a rational polynomial function of unknown parameters, or on assuming a specific time function input u(t) of unknown parameters. In simple control problems, (a second-order system with one or two controller parameters) one can apply direct computation to obtain an explicit optimal solution in terms of these parameters. On the other hand, one can apply the available finite-dimension optimization methods [14, 15] numerically to solve a general control problem. In this context, a technique called a control vector parameterization is on hand. The second approach is to derive a control law based on certain strategy without assuming any prior controller structure or function. In other words, the optimal controller structure or function and its parameters are determined in one solution framework. This infinitedimensional optimization is the kernel of the optimal control subject.

1.4. FINITE DIMENSIONAL OPTIMIZATION–A MATHEMATICAL REVIEW A brief review of finite-dimensional static optimization will be given in this section. It is useful to have some knowledge, as there are similarities, and as we borrow often some of the developed results of the finite-dimensional optimization. From such a rich subject, we shall introduce only the most important concepts, definitions, and theoretical results assuming that the reader knows the basic principles of unconstrained minimization or maximization of single variable scalar functions. Typical objective is to minimize a nonlinear function f(x) (objective function) or to maximize – f(x) where the variable x is in general multidimensional [x1 x2 x3… xn] vector. The point of extreme is denoted as x*. In such a framework, we use the terminology of nonlinear programming NP. The NP is applied to both of unconstrained and constrained functions. Let us first define the two concepts of minima: the strong minimum and the weak minimum, and the concept of convex functions.

Mathematical Background and Optimal Problem Modeling

5

Definition 1.1 - strong minimum A point x* is said to be a strong minimum of the function f(x) if a scalar α > 0 exists such that

f ( x*)  f ( x * x)  x  0 || x || 

(1.1)

It simply means that the objective function increases locally in all directions. Definition 1.2 - weak minimum A point x* is said to be a weak minimum of the function f(x) if a scalar α > 0 exists such that

f ( x*)  f ( x * x)  x  0 || x || 

(1.2)

It means that the objective function remains unchanged in some directions and increases locally in other directions. In reference to these definitions, we said that the function possesses a unique global minimum if these definitions hold for α = ∞ or otherwise a local minimum. Definition 1.3 - convex and strictly convex functions If f(x) is a scalar valued function of a vector x, it is said to be convex if, for distinct vectors x1 and x2 and a real scalar α, the following relation holds:

 f ( x1 )  (1   ) f ( x2 )  f ( x1  (1   ) x2 ), 0    1

(1.3)

A strictly convex function has only the inequality sign in Equation (1.3). Strictly convex function has at most, one minimum. However, not every strictly convex function has a minimum. Concerning the unconstrained case, the derivation of the necessary and sufficient conditions for a minimum is the most important theoretical results. For continuously differentiable function f(x), there are two necessary conditions derived from the approximation of f(x) by Taylor‟s expansion around an arbitrary neighborhood point, i.e. n

f ( x   )  f ( x)    i i 1

f 1 n n 2 f    i  j | x  ... xi 2 i 1 j 1 xi x j

(1.4)

where 0< ζ < 1. If, moreover, f(x) is a strictly convex function then - First-order necessary condition of optimality is



f ( x*)  0, f  f x

1

fx

2





f xn ,

fx 

f x

(1.5a)

6

Optimal Control Engineering with MATLAB

The symbol  stands for the gradient and the point x* satisfying this condition is called a stationary point. A second-order necessary condition for the strong minimum is

 f x1x1 f x x 2 2  f ( x*)  0,  f :  2 1     f xn x1

fx x



fx



1 2 2 x2

 f xn x

f x xn  1 f x xn  2 ,    f xn xn  

 2



f xy 

2 f xy

(1.5b)

In other words, the matrix  2 f should be positive definite. This condition must be used in combination with the first order necessary condition. In addition, we can distinguish the minimum case from the maximum one by the inequality sign. - A second-order sufficient condition for  2 f  0 , so x* is actually a strict minimum. Optimization problems can be either with equality or inequality constraints. Let us consider first the problem of minimizing the function f(x, y) subjected to the equality constraints (In control problems, we usually have the variable x represents a system state vector and y represents either a manipulated variable vector or an input control vector u).

y   y1

g i ( x, y )  0,

y 2  y m , i  1,2,..., m, m  n

(1.6)

In special cases when xi states are available explicitly from the constraints in terms of y, then we can eliminate them from f and proceed as in the unconstrained case. However, in general we minimize instead a new scalar function as

H ( x, y,  )  f ( x, y )  T g ( x, y )

(1.7)

where the vector λ is an arbitrary constant (the Lagrangian multiplier). It can be shown that the ith Lagrangian multiplier is the negative sensitivity coefficient of the minimum value of f to a change in the ith constraint. The necessary stationary conditions of are

H  0, xi

H  0, i  1,2,..., n; g  1,2,...m y j

(1.8)

The positive definite condition (sufficient condition) is given by

H yy  H yx f x

1

f y  f y ( f x ) 1 H xy  f y ( f x ) 1 H xx f x f y  0 T

T

T

T

1

(1.9)

Optimization problems, including inequality constraints, occur frequently in control engineering, where upper and lower bounds on variables should be considered. For example,

Mathematical Background and Optimal Problem Modeling

7

the armature current or the shaft speed of a DC servomotor should not exceed some threshold values or the temperature must be within specified values, etc. Constraints can be mathematically represented by

g i ( x, y )  bi ,

y   y1

 y m , i  1,2,..., m, m  n

y2

(1.10)

Basically, by introducing a vector of real values (usually they are called slack variables) ν = [ν1 ν2 … νm ]T such that

 i 2  bi  gi ( x, y )

(1.11)

Then it is possible to redefine the optimization problem as follows: Minimize : f ( x, y)

~( x, y )  g ( x, y )    b, Subjected to : g 2

b  b1 b2  bm 

T

The so-called Kuhn-Tucker conditions [16] constitute the necessary stationary conditions for a minimum. It states

f g  T  0, xi xi

i  1,2,  , n

f g  T  0, y j y j

j  1,2,. , m

 j (b j  g j )  0,

j  1,2,  , m

 j  0,

j  1,2,  , m

g j  bj ,

(1.12)

j  1,2,  , m

Although the problem was formally converted into an equivalent case of equality constraints, the solution requires several checking possibilities. It is proven that if the function f(x, y) is strictly convex and the constraints g(x, y) are also convex then the Kuhn-Tucker‟s conditions are necessary and sufficient for a unique minimum of the function subject to the inequality constraints. The above approach solves the problem indirectly by utilizing the stationary conditions. Alternatively, the problem can be solved directly by employing different numerical iterative methods. The most elementary method is called steepest descent method, which uses the firstorder gradients. For minimizing unconstrained n-variables function f(x), the iterative procedure, which constitutes the method of the steepest descent, is

xi  xi j

j 1

k

f xi

j 1

, i  1,2,..., n ;

j  1,2,3,....

(1.13)

8

Optimal Control Engineering with MATLAB

where the initial n-component variable is x0 =[ x10 x20 …xn0]. The value k is an arbitrary constant, which is chosen (for all i) to give a suitable small step. The value of the step k should be carefully chosen. A large step may overshoot the minimum, whereas a small value may slow the convergence. However, as the minimum is approached the derivative approaches zero and hence the convergence becomes even slower, which is a great drawback. In fact, practically the steepest descent method is of little use due to this drawback. To avoid a slow convergence in the vicinity of a minimum and/or to reduce the computation burden and memory requirements, many other direct methods exist. For instance, method uses second-order gradient or conjugate direction method, which is developed for quadratic functions, Fletcher-Reeves method, which is developed for smaller memory requirements, Powell's method having the property of quadratic convergence and does not calculate the gradient, etc. [17]. In optimization MATLAB toolbox, the fminsearch is a multidimensional unconstrained nonlinear minimization function. This function uses a direct search method called simplex search, which does not use numerical or analytical gradient [18]. However, this function may only give local solutions over a real domain; for complex parameter vector, they must be split into real and imaginary parts. In chapter two, the fminsearch function will be used to solve the parameter optimization problems in control systems. Direct methods for equality-constrained minimization were also developed. For the equality-constrained case, there are mainly two approaches. The first approach is to determine the x vector by solving numerically (or analytically) the constrained Equations (1.6) for given vector y and then to implement one of the optimization methods mentioned above, i.e. the iterative procedure starts with initial y vector to determine the initial x vector and consequently, starts the optimization method. The second approach is to lump the given function and the constraints in one function say H (Equation (1.7)) using Lagrangian multipliers, λ, vector given by

T  (

f T g 1 ) [ ] x x

(1.14)

The inequality constraints are further complicating the search for an optimum. One efficient technique is to use the transformation of the inequality constraints to the equality constraints. For example, some constraints and their corresponding equality transformations [19] are

1  y  1  m1  y  m2

y  (2 /  ) tan 1 ( x)

 y  m1  (m2  m1 ) sin 2 ( x) y  0  e x or x 2

(1.15)

However, in some cases, such transformations may complicate the problem or the computation burden.

Mathematical Background and Optimal Problem Modeling

9

1.5. OPTIMAL PROBLEM MODELING The first step in controlling physical systems (plants) is to have a mathematical model that describes the system dynamic as accurate as possible, and in the same time the model should not be so complicated as well. A mathematical model is an aggregation of equations and numerical data, which are employed to describe the dynamic behavior of a system in quantitative terms. Among several approaches of mathematical modeling, we are interested in the state-space models as it is the most suitable approach for optimal controller design; differential equation models or transfer function models can be easily transformed to statespace models. In this book, the dynamic systems are those governed by differential (continuous-time systems) or deference (discrete-time systems) equations. Any set of differential equations can be put into the non-autonomous or the time-varying form

dx  x(t )  f ( x(t ), u (t ), t ) dt

(1.16)

where the function f is a linear or nonlinear function, and the state (state variables) vector and the control (control variables) vector are

x(t )  x1 (t ) x2 (t )   xn (t )   n , T

u (t )  u1 (t ) u 2 (t )   u m (t )   m , m  n T

If the time does not explicitly appear in the function f, then the so-called autonomous or time- invariant model is obtained. When the plant is composed of subsystems of different order differential equations then conventionally we convert to a first-order differential set. This can be performed by converting the n-order differential equation to n first-order differential equations using a controllable (or observable) canonical form. For example, a plant consisting of one 1st-order, two 2nd-order, and one 3rd-order subsystems can be written by the form 1.16 with a dimension of 8. The state-space model should relate the physical quantities and their derivatives to the input by a minimum number of states, which exactly describe the system behavior. On the other hand, the state set is not unique and the selection of specific set depends on the problem formulation. One of the biggest advantages of these models is the mathematical generalization. This means that different physical systems can have the same state-space model. Although we have not yet defined precisely the optimal control problems, it is advantageous to model some physical systems in state-space terminology and to suggest specific optimal control problems. Consequently, the given examples demonstrate the practical aspect of the optimal control theories that will be considered later.

10

Optimal Control Engineering with MATLAB

Example 1.1. A pure double integral system It is the simplest second-order system. The model is often used for preliminary analysis of many complicated practical systems such as the vehicles, rockets, the drive of the hard disc, some special types of AC motors, etc. Such a model may represent a simple transitional motion such as an object sliding on a frictionless surface. A rotational motion such as orientation of satellite with thrusters may be also simply modeled by this system. Let us model a car as a point mass m with the acceleration a(t) and driving (braking) force F(t) acting on it such to move a distance d(t) on a frictionless road as shown in Figure 1.1.

Figure 1.1. A mass point car model.

Newton's second law of motion states

F (t )  ma (t )  m d (t )

(1.17)

For this system, let the state variables be defined as the position and the velocity of the car, i.e. x1 (t )  d (t ), x2 (t )  d (t ) , and the ratio F(t) / m be defined as the control variable u(t). If we assume that the car was initially at rest, then the double-integral model is

x1 (t )  x2 (t ), x1 (0)  0 x2 (t )  u (t ), x2 (0)  0

(1.18)

One basic optimal control problem is: to determine the force function behavior which transfers the car in minimum time tf to a point a D distance away and where the car has S velocity, i.e. to reach the point defined by the terminal state variables x1(tf) = D and x2(tf) = S. In this case, the optimal criterion is tf

min{J (u )}  min{  1.dt }  min{ t f }

(1.19)

0

Another problem is to reach specific terminal states with a minimum time and minimum control effort (fuel) as well. In this case, the optimal criterion has the form

Mathematical Background and Optimal Problem Modeling

11

tf

min{J (u )}  min{  ( 1   2 u 2 ) dt }

(1.20)

0

where γ1 and γ2 are some appreciated weights. These weights should be chosen without exceeding the maximum possible available value of the used force. A third problem is to assume that the car has an initial speed, and it is required to stop the car after T units of time and at D distance using a minimum amount of fuel. In this case, we can assume without loss of generality that at t = 0, x1(0) = 0. The optimal control problem is as follows:

 x1 (t )  x 2 (t ), x1 (0)  0, x1 (T )  D  x 2 (t )  u (t ), x 2 (0)  S , x 2 (T )  0

(1.21)

T

min{ J (u )}  min{  u 2 dt } 0

There are some other problems, which can be also addressed including constrained problems. For instance, car velocity should not exceed a specified value or to satisfy certain time function. Furthermore, other various systems can have a similar state-space model, and hence similar optimal control problems can be defined. Intuitively, the typical plant specifies the number and types of optimal control problems. On the other hand, in general, the optimal control theory does not depend on the control system dynamic. In fact, the optimal control theory is equally applied to the time-invariant, time-varying, linear, and the nonlinear systems.

Example 1.2. Armature DC servomotor An armature DC servomotor that derives an inertia load is shown in Figure 1.2.

Figure 1.2. An armature DC servomotor.

Optimal Control Engineering with MATLAB

12

The symbols are defined as follows:      

Ra and La are the armature resistance and inductance respectively. ea and eb are the armature and back-emf voltages respectively. ia is the armature current. ζm , wm are the angular position and angular velocity at the load side respectively. τm and τl are the motor and load torques respectively. J and B are the load moment of inertia and viscous friction of bearings respectively.

The task first is to derive a mathematical description in state-space form. Assuming that the servomotor is at rest initially, the following four equations describe the relations between the physical quantities: - The armature electric circuit is described by

e a  R a i a  La -

dia  eb dt

The back-emf of the motor, which is assumed to be linearly proportional to angular velocity, thus

eb  K b w m -

(1.23)

The motor torque characteristic, which is assumed to be linearly proportional to armature current, thus

 m  K i ia -

(1.22)

(1.24)

The rotor dynamic of the load; the load is represented by a torque and bearing friction, thus

 m   l  B wm  J

dwm d , wm  m dt dt

(1.25)

where Kb and Ki are the two proportional constants. It can be noted that the variables, which change dynamically with time are, ia, wm, and ζm, and hence we should choose them as state variables. Therefore, we select the three states as

x1 (t )   m (t ) x 2 (t )  wm (t )

(1.26)

x3 (t )  ia (t ) By substituting in the above equations, we can write the following third-order state- space model:

Mathematical Background and Optimal Problem Modeling

13

x1 (t )  x 2 (t ) x 2  (t ) 

1 ( Bx 2 (t )  K i x3 (t )   i ) J 1 x3 (t )  ( K b x 2 (t )  Ra x3 (t )  ea ) La

(1.27)

Let us also define the following for the system input and output: - Two inputs, which are u(t) = ea and w(t) = τl, which is considered as a constant deterministic disturbance. -

One output y(t) = ζm(t)

Therefore, we can put all these equations in the following matrix model:

  x   0 1  1   B  x 2   0  J    Kb  x3   0  L a 

     0  0   x1   0  Ki      1 x  0 u  2  J  w J    1   0  R   x3      L  a  a  La 

y   m  1 0 0 X ; X  x1

x2

(1.28)

x3 

T

Conventionally, we add a dummy state to represent the constant deterministic disturbance input (or even input of ramp function). The constant load torque is represented by a fourth state x4, and in turn this fourth state has zero derivatives (if a ramp disturbance is considered then two dummy states should be added). The state-space model in 1.28 becomes

1  x   0 B 1    0     x2  J     Kb  x 3  0  La    x   4  0 0

0 Ki J R  a La 0

y   m  1 0 0 0 X ;

0  0 1   x1        0  J  x2      1  u 0   x3   L   x   a  0   4   0  X  x1

x2

x3

(1.29)

x4 

T

We assume that the control of the servomotor should be performed within the time interval [0, tf]. It will be assumed that the final time is much greater than the largest motor time constant, i.e. a steady-state case. It is to find a state feedback control law in the form u = KX, where K is a constant gain vector, such that minimum energy dissipation in the motor resistance Ra is secured, and a minimum supplied input energy is obtained during a specific or

14

Optimal Control Engineering with MATLAB

accepted regulating time (rise time and/or settling time). Therefore, a performance index describing these demands is given by 

J   ( Ra ia  rea ) dt 2

2

(1.30)

0

The weighting constant r is adopted to adjust the required transient and steady- state responses (or to reach as possible acceptable responses). Alternatively, in terms of model states as well as the system input one can put the performance index to the form

0   T 0 J   (X 0 0  0

0 0 0 0 0 Ra 0

0

0 0 X  ru 2 ) dt 0  0

(1.31)

Of course, u represents the output of the optimal controller or in other words, the reshaping of the armature voltage ea. A certain value can be assigned to the (4, 4) element within the matrix in J to have the possibility of accounting for the load torque effects. This problem as stated in Equations (1.30) and (1.31) is one of the typical optimal problems so-called linear-quadratic optimal control, see chapter seven. Such problems include also the case where the final time tf is somehow comparable to the system time constants. In some other applications, it is required to maintain the motor speed constant at a particular value by controlling the input voltage. In this case, the system becomes of a second-order. In a similar way, models of DC shunt or series and AC motors, which are usually found in industrial applications, can be formulated.

Example 1.3. Double-mass spring-damper system In many mechanical engineering applications, the designer must carefully study the structural resonances that usually occur. The double-mass spring-damper system shown in Figure 1.3 can be used for such studies. The spring and damper coefficients of this mechanical compliance are k and b respectively. The system input, u(t), is the external force f (the controller output), while the displacement y2 is the system output to be controlled. It is assumed that the spring satisfies Hook‟s law. The damper is linear to the velocity, and the sliding surface is frictionless. Applying the Newton‟s second-law, the equations describing the movement of this system can be written as

Mathematical Background and Optimal Problem Modeling

   m1 y1  f  k ( y1  y 2 )  b( y1  y 2 )    m2 y 2  k ( y 2  y1 )  b( y 2  y1 )

15

(1.32)

In this problem, we have two second-order differential equations, which mean that we need to define four state variables to describe this system; two positions and two velocities of the masses, m1 and m2. We select a vector of four states as

 x1 (t )   y1 (t )   x (t )  y  (t )  X (t )   2    1   x3 (t )   y 2 (t )        x 4 (t )  y 2 (t )

(1.33)

Figure 1.3. A double-mass spring-damper system.

Then a state-space model in matrix form is given by

1 x    0 k b 1        x2  m1 m1     0 0  x3   k b x     4   m2 m2 y  0 0 1 0 X

0 k m1 0 k  m2

0  b   x1   10      m1   x 2   m   u  1   x3   01  b      x   m2   4   0 

(1.34)

For specific masses, the spring and damper coefficients constitute the two factors affecting the system natural frequency and damping ratio. In most practical applications, even with simple feedback controller the open-loop unit-step response is unstable. Moreover, the values of k and b change due to many reasons such as temperature or liquid properties, etc. Therefore, a robust controller is essential. Any designed controller should stabilize the system and satisfy fast settling time and minimum overshoot of the output (the displacement y2). The

16

Optimal Control Engineering with MATLAB

design method so-called pole-placement can be applied to obtain the required controller, but a sort of trial-and-error procedure should be adopted to reach precise values of specifications. In some applications where the minimum terminal time, minimum input energy, and minimum terminal values of error are demanded, the optimal criterion becomes tf

min{ J (u )}  min{[ X (t f )  X r ]T M [ X (t f )  X r ]   (1  ru 2 ) dt }

(1.35)

0

where M is a positive definite symmetric matrix, X(tf) and Xr are the terminal and required state vectors respectively, and r is the energy efficiency factor.

Example 1.4. Van der Pol oscillator In optimal control studies, the Van der Pol oscillator is commonly used to illustrate various optimal problems including time-optimal control, singular control, and linearquadratic optimal control. Moreover, since the oscillator is a nonlinear dynamic system various numerical methods have been applied to study their characteristics. Figure 1.4 depicts a modern RLC semiconductor model.

Figure 1.4. The semiconductor Van der Pol oscillator.

A differential equation describing the oscillator circuit is given by a 2nd-order nonlinear differential equation of the form

v0   (  v0 )v0  v0  vi ,   0 2

(1.36)

where μ and ε are the parameters that controlling the shape of oscillation. The values of these parameters are determined by the values of the passive elements, while the diode characteristics specify the nonlinear term. Let us define the states and the input as

Mathematical Background and Optimal Problem Modeling

17

x1 (t )  v0 (t ) x2 (t )  v0 (t ) u (t )  vi Therefore, a state-space model is

 x1  x2 , x1 (0)  x10  x2   x1   (  x12 ) x2  u, x2 (0)  x20

(1.37)

For given initial conditions, a basic optimal problem is to require a minimum time control, and at the same time satisfying terminal constraint as well as certain input constraint, specifically tf

min{  1.dt} subjected to 0

(1.38)

x12 (t f )  x 22 (t f )  m 2 ,  M  u (t )  M , t  [0, t f ] The optimal solution is a bang-bang control (either u = - M or u = + M), which responds rapidly and satisfies the terminal condition. On the other hand, if, for example, we cancel the terminal constraint and replace the performance index by a quadratic one and for specified value of terminal time, i.e., tf

J (u )   ( x12  x 22 ) dt, t f  T ,

(1.39)

0

then we will have a special optimal control problem, which is called a singular control problem. Furthermore, another non-singular control problem is solved in Example 4.9. Unconstrained problems with finite or infinite terminal time can also be stated with the quadratic performance index to derive a numerical optimal solution u(t) trajectory or a u(x1,x2) state feedback control law.

Example 1.5. Two-link robot manipulator (2-Link arm) This is another important example that is used to demonstrate the derivation of many control strategies, including optimal control, adaptive control, robust control, etc. Both unconstrained and constrained cases can be considered when an optimal control law is required. Due to the relatively higher order and system non-linearity, a closed-form solution is impossible, and hence it is a very good example to test the characteristics of many numerical techniques.

Optimal Control Engineering with MATLAB

18

There are several approaches to derive a mathematical model of robot manipulators. However, a usual approach is to use Lagrange's equation of motion, which relates generalize forces to generalize coordinates through the potential P and kinetic K energies of a conservative system. Let us consider only the motion of the manipulator in its vertical plane. We denote the difference (K – P) by the Lagrangian function L.

LK P

1 2 mv  mgh 2

(1.40)

where m is the mass, v is the velocity, h is the height, and g is the gravitational field acceleration. Lagrange's equation states that

d L L ( )  , dt q  q

(1.41)

where q is an n-dimensional vector of generalized coordinates qi , and τ is an n-dimensional vector of generalized forces τi.

Figure 1.5. Two arms robot manipulator.

Let the first-quadrant represent the plane of the workspace of the two links manipulator as shown in Figure 1.5, where m1, m2, l1, and l2 are the masses and lengths of the two arms respectively. The joint angles are ζ1 and ζ2 (ζ1 for shoulder, ζ2 for elbow).We shall assume that the masses of the arms are uniformly distributed over the length, the arm moments of inertia, end-effector mass and the joint friction is negligible. For such a manipulator, it can be shown that the Lagrangian function is

L  K1  K 2  P1  P2 L

1 1 m1l121 2  m2 (l121 2  l 22 (1   2 ) 2  2l1l 2 (1   2 ) cos 2 ) 2 2  g (m1l1 sin 1  m2 l1 sin 1  m2 l 2 sin(1   2 ))

(1.42)

Mathematical Background and Optimal Problem Modeling

19

Carrying out the required partial and ordinary differentiations in Equation (1.41) yields the following matrix equation so-called Euler-Lagrange equation:

       M (1 , 2 )  1   C (1, 2 ,1 , 2 )  1   G(1 , 2 )   1   2   2   2 

(1.43)

where the matrices M, C, and G are given by

m l 2  m2 (l12  l 22  2l1l 2 cos 2 ) m2 (l 22  l1l 2 cos 2 ) M  11  m2 (l 22  l1l 2 cos 2 ) m2 l 22  

   C  m2 l1l 2 sin 2  2  1

 1   2   0 

(m  m2 )l1 cos1  m2 l 2 cos(1   2 ) G  g 1  m2 l 2 cos(1   2 )   Equation (1.43) can be rewritten as

  1 1  1  1 1  1     M C     M G  M    2  2  2 im11 im12  1 where we denote the matrix M    im21 im22 

(1.44)

A vector of four states should be utilized to describe the dynamic of the robot manipulator. We select the state vector, which is composed of the joint angles and velocities, i.e.

x1

x2

x3

  T x 4   1  2 1  2   

T

In terms of the selected states, the elements of the matrix M-1 can be determined. They are given as

m2 l 22 m (l 2  l l cos(x 4 )) , im12   2 2 1 2 dm dm 2 2 m (l  l l cos(x 4 )) m l  m2 (l12  l 22  2l1l 2 cos(x 4 )) im21   2 2 1 2 , im22  1 1 dm dm im11 

Optimal Control Engineering with MATLAB

20 where,

dm  det(M )  m 2 l12 l 22 (m1  m2 )  2m22 l1l 23 (cos(x 4 )  1)  m22 l12 l 22 cos2 ( x 4 )

c   m l l x sin( x4 )  m2 l1l 2 sin( x4 )( x1  x2 ) c C   11 12    2 1 2 2  0 c21 c22   m2 l1l 2 x1 sin( x4 )   g   g ((m1  m2 )l1 cos(x3 )  m2 l 2 cos(x3  x4 )) G   11     gm2 l 2 cos(x3  x4 )  g 21    Now a state-space model can be written as

 x1   w1  x   w  2   3  x3   1     x4   0

w2 w4 0 1

0 0 0 0

0  x1  im11 im12  0  x2  im21 im22   u1   0  x3   0 0  u2      0  x4   0 0 

(1.45)

where,

w1  (im11c11  im12 c21 ), w2  (im11c12  im12 c22 ) w3  (im21c11  im22 c21 ), w4  (im21c12  im22c22 )

u1   1  g11 , u2   2  g 21 The above model does not include the dynamic of the coupled actuators. However, it can be easily modified to include the dynamic of the actuators; for more elaboration, the reader is advised to see the reference [20]. In addition, clearly, although the robot manipulator dynamic model has the property of linearity in parameters it is a coupled nonlinear multivariable system. Therefore, the solution of the dynamic model can be carried out only numerically, provided that the matrix M-1 exists at each step of integration. Basically, there are two main motion tasks of any robot manipulator in its workspace, the point-to-point, and the tracking motions. In the former task, the end-effector of the robot manipulator moves from a certain initial position (usually in Cartesian coordinates) to another position, which consequently becomes an initial position for next motion. The second motion describes the continuous tracking of the end-effector to an arbitrary trajectory (usually defined by certain dynamic equation). Accordingly, different performance indices may be found in optimal control. A minimum terminal time may be required with a point-to-point motion with input constraints, i.e. to solve the dynamic system for specific parameters of the used actuators and sensors. On the other hand, the tracking motion requires the minimization of the performance index.

Mathematical Background and Optimal Problem Modeling

21



J   (e(t ) T Q e(t )  u (t ) T R u (t )) dt

(1.46a)

0

The error e(t) is a vector defined by the difference between the actual angles and velocities of the joints and the angles and velocities of the required trajectory. The weights Q and R are suitable positive semi-definite and positive definite matrices. In 2-link motor control problems [21], the performance index is given as T

1 1 J  [ (T )   *)]T M [ (T )   *]   u (t ) T Ru (t )) dt 2 20

(1.46b)

where ζ is the joint angle vector and ζ* is the desired final posture. The first term in the index 1.46b means that the joint angle is going to the target value ζ*, which represents the reaching movement. The second term defines the energy efficiency.

Example 1.6. Lateral aircraft autostabilizer The aircraft motion is usually described by two mathematical models, the longitudinal, and the lateral. Both models are complicated, nonlinear, and multivariable dynamic systems. Moreover, based on the various flight modes and the assumed aeronautic conditions different models are obtained. We shall not go into such model elaboration, but instead we focus on the lateral motion that assumes a level flight path, a still air flight path, and a rigid body aircraft, see Figure 1.6.

Figure 1.6. Aircraft coordinates system.

Optimal Control Engineering with MATLAB

22

Let us consider the lateral motion defined by the following six state variables:  x1 is the sideslip velocity β  x2 is the roll or bank angle φ  x3 is the roll rate p  x4 is the yaw rate r  x5 is the aileron angle δa  x6 is the rudder angle δr For primary design, a linearize state-space model of the lateral motion is usually adopted and similarly, for their electro-hydraulic actuators. Such a model is given by

 x1   Yv  x   0  2   x3   L    x 4   N   x5   0     x 6   0

g /U 0 Lr

0 1 Lp

 1 Y a 0 0 0 L a

Nr 0

Np 0

0 0

N a  k1

0

0

0

0

  x1   0  x   0  2    x3   0    N  r   x4   0 0   x5  k 3     k 2   x 6   0 Y r 0 L r

0 0  0   u1   0  u 2  0  k 4 

(1.47)

where the airplane parameters 

Yv , Yδa , Yδr are the side force term derivatives due to sideslip, aileron and rudder angles respectively.  Lβ , Lr , Lp , Lδa , Lδr are the rolling moments due to the sideslip, yaw rate, roll rate, aileron and rudder respectively.  Nβ , Nr , Np , Nδa , Nδr are the yawing moments due to the sideslip, yaw rate, roll rate, aileron and rudder respectively.  u1 , u2 are two forcing electrical inputs supplied to the electro-hydraulic control surface actuators.  ki , are the actuators parameters. In this example, additional state variables must be introduced to correspond to a reference variable (as in Example 1.2 when one additional state is introduced to correspond to the load disturbance). The reference variable could be a demanded roll-rate or yaw-rate. Conventionally, the design of the auto-stabilizer controller demands a specific response, which does not necessarily follow the reference variable). In fact, the response is usually assumed to be that one which corresponding to a second-order dynamic defined by specific damping ratio and natural frequency. Therefore, if a short impulsive command for roll-rate is applied, the approximate response will be described by the state equation

x7  x8 x8  a1 x7  a 2 x8 , a1 , a 2  0 where x7 is the approximate roll-rate, and a1, a2 are constant design parameters.

(1.48)

Mathematical Background and Optimal Problem Modeling

23

Let us assume that the roll-rate command is applied at t0, and we require stabilizing the aircraft laterally after specified units of time tf such to:  Minimize the sideslip velocity x1.  The actual roll-rate x3 is approximating the desired response x7, i.e. minimum error.  Minimize the control effort, i.e. restrict the actuator signals u1 and u2. An optimal criterion, which places emphasis on the above, is given as t0  t f

min{ J (u1 , u 2 )}  min{

t (

1

x12  ( x3  x7 ) 2   2 u12   3 u 22 ) dt }

(1.49)

0

where γ1 , γ2 , and γ3 are appreciated weights.

Example 1.7. An inverted pendulum mounted on a cart This is another example, which is usually used in modern control theory. It is invoked to elucidate the application of different techniques, including optimal control, adaptive control, robust control, etc., [22]. The mathematical model can be derived on the same theoretical bases that be used for robot manipulator. Figure 1.7 shows an inverted pendulum mounted on a cart, which is moved on a frictionless surface, where:  The mass of the cart is M.  The mass of the pendulum is m.  The length of the pendulum is L.  The displacement of the cart is y.  The angular position of the pendulum is ζ.  The gravitation is g. The force acting on the cart (by the motor inside the cart) is u, which is the input to the system.

Figure 1.7. Inverted pendulum mounted on a cart.

Optimal Control Engineering with MATLAB

24

Similar to what was done in Example 1.5, the Lagrangian function is given by

1 1 L  K  P  ( My2  m 2 )  mgl cos( ) 2 2

(1.50)

where y ,   are respectively the velocity of the cart and the angular position of the pendulum. The pendulum velocity has two components, a vertical vz and a horizontal vy, thus one can write the following dynamic:

d d [ y  l sin( )) 2  ( l cos( )) 2 dt dt 2 2  ( y   l  cos( ))  (l  sin( ))

w 2  v y  v z  ( 2

2

(1.51)

 y  2  l 2  2  2lx   cos( ) Two Lagrangian equations for the cart and the pendulum motions should be considered.

(1.52)

d L L ( )  Fy dt y  y

(1.53)

Working out the derivatives and substituting in the Equations (1.52), and (1.53) yields

y  

ml  2 sin( )  mg sin( ) cos( )  u M  m sin 2 ( ) (1.54)

  

 ml  2 sin( ) cos( )  ( M  m) g sin( )  u cos( ) l ( M  m sin 2 ( ))

If a state vector is defined as

x1 Then the state equations are

x2

x3

x4    y T

y   

T

Mathematical Background and Optimal Problem Modeling

25

x1  x2 mlx42 sin( x3 )  mg sin( x3 ) cos(x3 )  u x2  M  m sin 2 ( x3 ) x3  x4 x4 

(1.55)

 mlx42 sin( x3 ) cos(x3 )  ( M  m) g sin( x3 )  u cos(x3 ) l ( M  m sin 2 ( x3 ))

The above nonlinear model can be linearized around a specific operating point, for example,   0 eventually    (the vertical-up position eventually, the vertical- down position). A small perturbation is usually assumed to initiate the pendulum motion, see problem 7.13. The standard optimal control problem is to determine the control that keeps ζ as close to zero as possible (minimal swinging) with minimal supplied energy and after elapsing finite time. Therefore, the performance index should be of the form tf

J (u )   (qx32  ru 2 ) dt

(1.56)

0

where tf is either a given specific value (could be much larger than the system time constant) or free, q and r are the relative weights. A large value of q corresponds to the angle ζ value, which is very near to zero, whereas a large value of r corresponds to energy saving. Another desire could be to return the disturbed pendulum arm to the vertical position and keep it unmoved there in minimum time and minimal supplied energy while the cart moves freely. In this case, different initial, terminal conditions and performance index have to be adopted, as given below.

x(0)  x10

x 20

x30



x 40 , x(t f )  x1 f

x2 f

0 0

tf

J   (   ru 2 )dt;   0

 (1.57)

0

where x1 f, x2 f are arbitrary cart terminal position and velocity respectively.

Example 1.8. A moon lander The task is to bring a spacecraft to a soft landing on the lunar surface using the least amount of fuel (maximizing the remaining fuel), Figure 1.8. Second Newton‟s law represents a simple dynamic model

mh    gm  a

(1.58)

26

Optimal Control Engineering with MATLAB

Figure 1.8. A spacecraft landing on the moon.

where h(t) is the spacecraft instantaneous height, v(t) is the velocity or height rate, m(t) is the mass of spacecraft, which is changing as fuel is used up, g is the gravity acceleration and a(t) is the instantaneous thrust assumed to be bounded as,

0  a(t )  1

(1.59)

Let the states x1(t), x2(t), and x3(t) represent the height, velocity, and mass respectively. The control input u(t) is the acceleration a(t). Let also a positive constant k that defines the fuel consumption coefficients. A nonlinear state-space model is

 x1 (t )  x 2 (t ) u (t )  x 2 (t )   g  x3 (t )  x3 (t )  ku(t )

(1.60)

Minimizing the amount of fuel used up means maximizing the amount of remaining once the vehicle has landed, and then an optimal criterion can be stated as

max{ m( )}  min{ m( )}

(1.61)

where τ denotes the first time (terminal time which is not given in advance) that both the height and the velocity are equal to zero, i.e., a terminal condition

x1 ( )  x2 ( )  0

(1.62)

Obviously, the optimal model should concern also the following two positiveness constraints:

x1 (t )  0, x3 (t )  0 It is an optimal problem of input and state constraints.

(1.63)

Mathematical Background and Optimal Problem Modeling

27

Example 1.9. A Thermal problem The problem is to cool a certain liquid in a tank (first) by pumping a cooling agent stored in an anther tank (second); Figure 1.9. Let us consider the following assumptions:      

The initial storage of the cold fluid is limited to C0 = 1. The instantaneous temperature of the first tank is x(t) whose initial value is equal to T0. The first tank has heat loss to the zero-temperature surrounding. The cold flow is described by a1 u(t), where u(t) is the control flow; a1 > 0. The heat carryover of the overflow liquid is a2 x(t) u(t); a2 > 0. The capacity of the transfer pump produces a constraint in the control flow u so that 0  u(t )  P .

Accordingly, the dynamic describing this thermal problem is simply given by a nonlinear first-order differential equation

x (t )   x(t )  [a1  a 2 x(t )] u (t ), x(0)  T0

(1.64)

One optimal problem may be stated as: obtain the optimal control u(t) that brings the temperature of the first tank down to zero in minimal time. tf

min {J (u )   1. dt}

0 u  P

(1.65)

0

Figure 1.9. Cooling a hot tank.

To make full use of the cooling agent, the optimal control can possess a coast-boost nature as follows:

Optimal Control Engineering with MATLAB

28

 0, 0  t  t1 u * (t )   P, t1  t  t f

(1.66)

where the times t1 and tf have to be determined. Therefore, the optimal control is to power-on the pump up to t1 units of time and then switched-off till the first tank temperature reaches zero after tf units of time; see problem 5.11.

Example 1.10. A container crane system A well-known control problem is to constrain the sway angle within a certain limit during load transportation in a container crane system. The crane system purpose is to lift and move loads through a hook suspended from a movable arm or trolley. Figure 1.10, shows an equivalent diagram for a suspended load on the electro-mechanical system, where      

Fy is the horizontal force acting on the trolley ML and MT are the load and trolley masses respectively yL and zL is the load coordinates in the vertical (y-z) plane ζ is the swing angle l is the length of the rope g is the gravitational constant

Let us assume a generalized coordinate q defined as q = [ζ, y]T, and a generalized input function U = [0, Fy]T. The kinetic energy T and the potential energy V functions are respectively given by

T

1 1 M T x2  M L ( y2  2 yl  cos  l 2 2 ) 2 2 V  M L g l cos (1  cos )

(1.67) (1.68)

The force Fy is generated by the DC motor whose input voltage is u. Similarly as in example 1.7, we can obtain the acceleration equations in the two general coordinates as follows:

y  

   

M L g sin cos  M L  2 sin  K M u ( M T  M L (1  cos2  ))

M L g sin  cos  M L  2 sin   K M u g cos  sin 2 l l ( M T  M L (1  cos  ))

where, KM is a motor constant.

(1.69)

(1.70)

Mathematical Background and Optimal Problem Modeling

29

Figure 1.10. The container crane system.

A linear model can be obtained by assuming a small swing angle such that the following approximations are valid: cos θ ≈ 1, sin θ ≈ θ, and θ2 ≈ 0. Therefore, we can rewrite equations (1.69, and 1.70) as

MLg K  M u MT MT

(1.71)

(M T  M L ) g K  M u MT l MT l

(1.72)

y 

   

A state-space model can be accommodated by assuming the following four states vector X as

X  x1

x2

x3

x4    y T

y   

T

In matrix form, the state-space model becomes

0  x1    x   0  2    x3  0   0  x4   

1

0 MLg 0 MT 0 1 (M T  M L ) g 0  MTl

0  0    x1   K M  0    x2 MT    u 0  x3   0  K 0  x4   M    M T l 

(1.73)

Usually, there are three essential practical aspects that describe the optimal performance. Safety and economic constraints require that both the load swing, and the transfer time are kept to minimum values during the load movement from an initial position to a desired specified position yd(tf). The following optimal criterion reflects these requirements:

30

Optimal Control Engineering with MATLAB

tf min{ J   t f   [ px22  q x42  r u 2 ] dt}

(1.74)

0

where α, p, q and r are appreciated positive constants.

All previous physical systems are represented by continuous dynamic models; however, there are many physical systems need to be modeled by discrete models. In fact, discrete models are essential for digital computer control. In chapter seven, some discrete issues are considered. The considered models represent a subset of a very large set, which includes many other application fields such as chemical, biological, financial, societal, etc. The idea was to model various kinds (electrical, mechanical, etc.) of systems for better understanding of the optimal problems. Nevertheless, the reader can later use the above examples to apply the optimal theories that will be presented in the next chapters. In many standard books of control theory, one can find different physical systems, which can be modeled as optimal control problems. In addition, some special control books consider the solution of a special optimal control problem; for example, the linear motion of electromagnetic levitators (LEL) system [23]. Finally, the reader can find many dynamic models in the comprehensive book of Klee [24], and hence various optimal control problems can be exercised.

PROBLEMS 1.1. Consider the problem of decelerating a car moving on a friction surface in such a way as to cover the least amount of distance in a given period T, with minimum energy loss. State an optimal problem, which assume an initial position of the car d, an initial velocity v0 and a friction coefficient b. 1.2. A boat travels with constant speed v and varying heading β through a region of strong currents. The motion from the shore is assumed to be in the x direction, while the current in the y direction with speed of w. It is required to cross the river and reach the opposite point in minimum time. Write a mathematical model of the boat motion and formulate the optimal problem. 1.3. Consider the suspension system is shown in figure 1.11. The springs satisfy Hook‟s law, and the damper is a linearly proportional to the velocity. The spring coefficients of the springs, the masses, and the coefficient of the viscous friction of the damper, are given constants. Our goal is to move the two masses to y1r and y2r respectively while minimizing a suitable quadratic cost of the position errors and energy for infinite control time. 1.4. A tunnel diode oscillator has an electrical circuit similar to that in figure 1.4. It is described by the Rayleigh equation (see problem 4.22), which is a second-order nonlinear differential equation of the circuit current and input voltage vi. In a fixed period T, it is

Mathematical Background and Optimal Problem Modeling

31

required to return the system to the equilibrium point, with minimum quadratic cost of the current and energy. 1.5 In a series RC-circuit derived by a source of V voltage, it is desired to charge the capacitance from zero initial value to a specific value q after T units of time, such that the power dissipation in the resistor is minimal. Model the optimal control problem. Hint: Assume that the state is the voltage across the capacitance. 1.6. Consider a cylindrical robot [25] with angular motion ζ and radial motion l’. The coordinated angular and radial motions in an assembly task, is shown in figure 1.12. The radial actuator force F and the angular actuator torque Tm are subjected to constraints. A mass m1 should be grasped by the robot at the initial position and moved to another position in minimal time.

Figure 1.11. A suspension system (Prob. 1.3).

Figure 1.12. A cylindrical robot (Prob. 1.6).

32

Optimal Control Engineering with MATLAB

1.7. A simple harmonic oscillator can be represented by the motion of a pendulum. Consider a pendulum in the vertical plane as shown in figure 1.13, where ζ(t) is the angular displacement. It is assumed that initially the pendulum has an angular displacement value p, and an angular velocity value w. Furthermore, assume that the angular displacement has small values, and the applied torque F(t) to the ball of the pendulum cannot exceed its mass in magnitude. The task is to bring the pendulum to rest as quickly as possible.

Figure 1.13. A pendulum in a vertical plane (Prob. 1.7).

REFERENCES [1]

[2]

[3]

[4]

[5]

E Junk, S. Lenhart, Z. Feng, "Optimal Control of Treatments in a Two-Strain uberculosis Model", Discrete and Continuous Dynamical System-Series B,Volume2, Number 4, 2002. Hem Raj Joshi," Optimal Control of an HIV Immunology Model", D.I.Cho, P.L. Abed and M. Parlar, "Optimal Production and Maintenance Decisions when a System Experiences Age-Dependent Deterioration", Optimal Control Appl. Meth. 14, 153-167 (1993). E. Esmailzadeh, H.D. Taghirad, "State-Feedback Control for Passenger Ride Dynamics", In the Transactions of the Canadian Society for Mechanical Engineering, 19(4):495-508,Dec 1995. Vesna Nevistic, James A. Primbs, "Constrained Nonlinear Optimal Control: A Converse HJB Approach", Technical Memorandum No. CIT CDS 96-021, December 1996. Chandeok Park, Daniel J. Scheeres, "Solution of the Optimal Control Problem using Hamiltonian Dynamics and Generating Functions", CDC Conference, Maui, Paper TuP07.3 (2003).

Mathematical Background and Optimal Problem Modeling [6]

[7]

[8]

[9]

[10] [11]

[12]

[13] [14]

[15] [16]

[17] [18] [19]

[20] [21] [22] [23]

33

Helmut Maurer, Inga Altrogge, Nadine Gorise, "Optimization Methods for Solving Bang-Bang Control Problems with State Constraints and the Verification of Sufficient Conditions", 44th IEEE Conference on Decision and Control, Seville, Spain December 2005. H. Maurer, D. Augustin, "Sensitivity Analysis and Real-Time Control of Parametric Control Problems using Boundary Value Methods ", On-line Optimization of Large Scale Systems, Spinger Verlag, Berlin 2001, pp. 17-55. T.C. Yang, H. N. Yu," Can you do it better? Design a Optimal Controller", R. Pfiffner, L. Guzzella, C.H. Onder," Fuel-Optimal of CVT Power Trains", Control Engineering Practice 11 (2003) 329-336. Feng Lin, Robert D. Brandt, "An Optimal Control Approach to Robust Control of Robot Manipulator", IEEE Transaction on Robotics and Automation, Vol.14, No. 1 February 1998. Mike Grimble, “Robust Industrial Control”, Prentice Hall Canada, 1997. Xuping Xu, Panos J. Antstklis, "Optimal Control of Switched Systems Based on Parameterization of the Switching Instants", IEEE Trans. Automatic Control, Vol. AC49, No. 1, January 2003. U. Ly, A. E. Bryson, and R. H. Cannon, "Design of Low-order Compensators using Parameter Optimization", In Applications of Nonlinear Programming to Optimization and Control, Laxenburg, Austria: IFAC 1983. P. M. Makila and H. T. Toivonen, "Computational Methods for Parametric LQ minimization: A survey", IEEE Trans. Automatic Control, Vol. AC-32, Aug 1987. Kuhn, H. W. and Tucker, A. W. , "Nonlinear Programming", Proceeding of the 2nd Berkely Symposium on Mathematical Statistics and Probability, Berkely, University of California Press, 1951. Singiresu S, Rao, "Engineering Optimization- Theory and Practice", Third Edition, New AGE International (P) Limited Publishers 1996. Lagarias, J.C., J. A. Reeds, M. H. Wright, and P. E. Wright, "Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions," SIAM Journal of Optimization, Vol. 9 Number 1, pp. 112-147, 1998. Box, M. J. "A Comparison of several current Optimization Methods and the use of Transformation in Constrained Problems", Computer Journal, 9, No.1, 1966. Mark W. Spong, Seth Hutchinson, M. Vidyasager, "Robot Modeling and Control", John Weily andSons, Inc 2006. Emanuel Todorov, Weiwei Li, "Optimal Control Methods Suitable for Biomechanical Systems", In 25th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, September 2003. Feng Lin, "Robust Control Design- An Optimal Control Approach", John Weily andSons, Ltd, 2007. Boldea I., Nasar S.A., "Linear Motion Electromagnetic Systems", New York, John Wiley, and Sons, Inc. 1985. Harold Klee, "Simulation of Dynamic Systems with MATLAB and Simulink", CRC Press 2007. Hans P. Geering, "Optimal Control with Engineering Application", Springer- Verlag Berlin Heidelberg 2007.

Chapter 2

CONTROLLER DESIGN BASED ON PARAMETER OPTIMIZATION 2.1. INTRODUCTION In this chapter, we shall discuss the design of fixed structure controllers, where the designer has to determine the unknown parameters of these structures. This approach is called the outward approach; the inward approach will be mentioned briefly. The controller parameters will be found according to a specific performance index. While tractable results are obtained for some indices, the big mathematical burden precludes the use of other indices. As with the design of the classical controllers, the transfer function concept is adopted to formulate the problem solution. The reader will note that the problem of finding the controller parameters is nothing than an ordinary optimization procedure. First, the concept of using the standard forms of the characteristic equation is highlighted. For this purpose, we use one of the famous standard forms, which are derived on the basis of the index of the integral-time-absolute error. For a thorough study, the applications to include the constrained case and the limitations of using this approach are discussed. Furthermore, we introduce another important approach named the control vector parameterization CVP. The CVP approach equally solves the problems of linear and nonlinear control systems. Finally, the genetic algorithm GA, which is one of the known evolutionary techniques, is used. Several examples are worked to explain the design procedure. The exact as well as the numerical solutions are considered. Moreover, to gain a view of comparison, the optimization procedure is applied to some fixed structure controllers, which were already designed by classical methods.

2.2. PERFORMANCE INDICES To meet a certain desirable set of performance specifications, one way to avoid the big burden of the trial-and-error procedure is to use single performance index IP to be maximized or minimized. The performance index is a number, which indicates the “goodness” of the system performance. Such IP could be established and based on which one may describe the goodness of the system response in a certain sense. In this way, the design procedure will be

Optimal Control Engineering with MATLAB

36

straightforward. Obviously, the use of the term “goodness” is somehow vague, unless the term “in a certain sense” is defined precisely. The goodness of the system response in a certain sense may not be so in another. Hence, the selection of certain IP requires itself an engineering skill, which deliberates for the meaning of the system variables. However, the selection of IP must offer good selectivity to distinguish non-optimal adjustment of the parameters, and it must exhibit a maximum or minimum value. In order to have a benefit tool, the IP must be a function of system parameters and must be easily computed. In control literatures, the usual used term is “optimal” which indicates goodness or best in a certain sense. The general form of these indices is an integration computed from initial time t0 (usually equal to zero) to final time tf (usually several times the longest time constant of the system). The integrand is a function of error signal and time, which appears explicitly or be implicitly in the integrand term. In the literature, various performance indices are proposed. The following indices are usually used: 



Integral square-error ISE  e (t ) dt 2

0





Integral-of-time-multiplied square-error ITSE  t e (t ) dt 2

0





Integral absolute-error IAE  | e(t ) | dt 0





Integral-of-time-multiplied absolute-error ITAE  t | e(t ) | dt 0





Integral square-error plus square error rate ISEER  [e (t )  e (t )] dt 2

2

0



Integral-of-time-multiplied absolute-error plus absolute error rate ITAEER 

 t [| e(t ) |  | e(t ) |] dt

0

There are also other indices such as the integral-of-square time-multiplied error. In practice, the choice of one of these indices depends on a tradeoff between the good selectivity (to obtain a precise value of the optimal parameters) and good insensitivity (a small variation of the index value over a range of values of the parameters). Obviously, the error signal can be replaced by any signal in the closed-loop control system; the concept of system state is used to define any signal. It is shown elsewhere [1] that for a second-order control system defined by a natural frequency ωn equal to 1 radian/second and a damping coefficient ζ , the above indices (denoted J) can be analytically computed in terms ζ. They have the following forms: 

For ISE index J   

1 4

,  0

Controller Design Based on Parameter Optimization

1



For ITSE index J   2 



For IAE index J  2



For ITAE index J  4  1



For ISEER index J   



For ITAEER index J  4  2  1 ,   1

37

,  0

8 2

,  1 ,  1

2

1 2

,  0

2

The plots of these indices with respect to the damping coefficient show that the ITAEER index has the best selectivity among the other indices, while the ISE index has relatively better insensitivity. It can be then easily computed that for ITAE index, the optimal value of the damping coefficient is 0.707. When δ = 0.7 is set, the second-order system results in swift response to a step input with approximately 5% overshoot. Although for high-order systems, these figures are different, the dominated second-order system can approximately express the transient and/or the steadystate responses. Since Graham and Lathrop introduced in their paper [2] the synthesis of optimum response, a great deal of interest in control system optimization is developed. For various order systems, they suggested the form of the closed-loop system transfer functions. These systems will achieve zero steady-state step (position) error systems, zero steady-state ramp (velocity) error systems or zero steady-state parabolic (acceleration) error systems and minimized ITAE. The results were found experimentally. For zero steady-state step error nth-order systems, the suggested closed-loop transfer function is given as

G(s) 

s  a1 s n

n 1

an  ...  a n 1 s  a n

(2.1)

Similarly, for zero steady-state ramp error systems, the suggested closed-loop transfer function is

G( s) 

a s  an s  a1s  ...  an1 s  a n n

n 1 n 1

(2.2)

Tables 2.1 and 2.2 show the optimum forms of the closed-loop transfer functions based on the ITAE index up to sixth-order for zero steady-state step and ramp input signals. These standard forms of transfer functions provide a simple and fast method for synthesizing an optimum dynamic response. Other forms based on other indices are also derived [3]. The closed-form results can be obtained for both of the following performance indices:

38

Optimal Control Engineering with MATLAB 

J 0   [ x(t )  x()]2 dt 0 

J 1   [ x 2 (t )  x  2 (t )] dt 0

where x(t) usually represents the system output (or the error actuating signal) and x(∞) is the corresponding steady-state value. The first index quantifies the speed of the control system for step input, and the second index is established to obtain best system tracking to exponential input of time constant τ. Let the output to a unit-step input is given by the rational transfer function

b0  b1 s  ...  bm s m 1 Y ( s)  a0  a1 s  ...  a n s n s

nm

(2.3)

where ai and bi are the system parameters. Then the first index J0 is computed and has the closed-form

J0 

1 ( B0  0  B1 1  ...  Bm  m  2b0 b1 ) 2a02 

(2.4)

where Δ is the determinant given as Table 2.1. Optimum Form of the Closed-Loop Transfer Function Based on the ITAE index for Zero Steady-State Step Error Systems

G(s) 

an an   nn n n 1 s  a1 s  ...  a n 1 s  a n

_________________________________________________________________________________

s  n _________________________________________________________________________________

s 2  1.4n s  n2 _________________________________________________________________________________

s 3  1.75n s 2  2.15n2 s  n3 _________________________________________________________________________________

s 4  2.1n s 3  3.4n2 s 2  2.7n3 s  n4 _________________________________________________________________________________

s 5  2.8n s 4  5n2 s 3  5.5n3 s 2  3.4n4 s  n5 _________________________________________________________________________________

s 6  3.25n s 5  6.6n2 s 4  8n3 s 3  7.45n4 s 2  3.95n5 s  n6 _________________________________________________________________________________

Controller Design Based on Parameter Optimization

39

Table 2.2. Optimum Form of the Closed-Loop Transfer Function Based on the ITAE index for Zero Steady-State Ramp Error Systems

G(s) 

a

s  a1 s n

s  an an   nn  ...  a n 1 s  a n

n 1 n 1

________________________________________________________________________

s 2  3.2n s  n2 ________________________________________________________________________

s 3  1.75n s 2  3.25n2 s  n3 ________________________________________________________________________

s 4  2.41n s 3  4.93n2 s 2  5.14n3 s  n4 ________________________________________________________________________

s 5  2.19n s 4  6.5n2 s 3  6.3n3 s 2  5.24n4 s  n5 ________________________________________________________________________

s 6  6.12n s 5  13.42n2 s 4  17.16n3 s 3  14.14n4 s 2  6.76n5 s  n6 ________________________________________________________________________

 a0  0   0   det(  0  ......   0

 a2 a1  a0

a4  a3 a2

 a6 a5  a4

... ... ...

0 ...... 0

 a1 ...... 0

a3 ...... 0

... ..... 0

    ) 0  .....   a n 1  0 0 0

(2.5)

The Δ j (j = 0, 1, 2… m) are determinants that are obtained from the determinant Δ by replacing the (j + 1)th column by the column [a1 a0 0 … 0]T. The coefficients B0… Bm are computed from the relations

B0  b02 B1  b12  2b0 b2 .......... .......... ........ Bk  bk2  2bk 1bk 1  ...  2(1) k b0 b2 k

(2.6)

.......... .......... .......... .......... .......... .......... ... Bm  bm2 The second index can be computed by two separate integrals. The first integral J0 is computed as above, while the integral J1 is computed by adopting one of the Laplace transform properties.

Optimal Control Engineering with MATLAB

40

y (t )  L1{sY ( s)} i.e. sY ( s) 

b0  b1 s  ...  bm s m a 0  a1 s  ...  a n s n

nm

(2.7)

When it is the case where e(t) = r(t) - y(t), then the above method (described by Equations (2.3-6)) yields exactly the same results as the ISE index gives. However, in control systems, the Laplace transforms of the error response for step input is usually given by n 1

i  bi s E ( s )  in0 i  ai s

(2.8)

i 0

Provided that the E(s) is stable then ISE index exists, i.e. it has a finite value. A closed form can be derived from the above-mentioned method. It is usually put in the following expression [4, 5]: 

2  e (t ) dt 

0

Hb , 2a n H a

(2.9)

where Ha is the determinant of the Hurwitz‟s matrix defined as

a n 1 a  n  0  H a  det( 0  0   0  

a n 3 an2 a n 1

a n 5 a n4 a n 3

.. .. .. .. .. ..

0 0 0

an ..

a n2 ..

.. .. a 0 .. .. a1

..

..

.. .. a 2

0 0  0  0 ) 0  a0   

(2.10a)

The determinant Hb is determined from the Hurwitz‟s matrix by replacing the first row by a row of elements, which are functions of the numerator coefficients of the error response (2.9) as follows:

[bn21

 (bn22  2bn1bn3 ) (bn23  2bn2 bn4  2bn1bn5 ) ... ... (1) n3 (b22  2b1b3  2b0 b4 ) (1) n2 (b12  2b0 b2 ) (1) n1 b02 ]

(2.10b)

The optimized parameters of the control system can be then obtained by the extreme of the closed form of the selected index in usual mathematical ways. The optimization of a single parameter is almost easy to carry out. However, when there are many parameters, this

Controller Design Based on Parameter Optimization

41

closed-form calculation becomes more tedious or may even results in unreliable parameter values. In such a case, the designer is postulated to assume non-optimal values of some parameters or to change the assumed structure of the controller and re-optimize the problem. This is the fundamental drawback of this approach.

2.3. PRELIMINARY DESIGN CONCEPT – THE INWARD APPROACH First let us consider a simple straightforward procedure design of a single-input singleoutput linear time-invariant control system. The design will be based on the optimum closedloop transfer functions; say, for example, that one based on optimizing ITAE index. Essentially, we will consider a control system configuration (plant-controller connection) and an optimum closed-loop transfer function of a certain order (inward approach). Concerning the classical feedback arrangement shown in Figure 2.1, a controller of a transfer function GC(s) is connected to the plant G0(s) in cascade and a unity-negative feedback completes the system configuration. It is the simplest configuration for utilizing the idea of inward approach. Let the required optimum closed-loop transfer function be denoted by T(s), which is given as

T ( s) 

GC ( s )G0 ( s ) 1  CC ( s )G0 ( s )

(2.11a)

Solving for the controller, we obtain

GC ( s) 

1 T ( s) G0 ( s ) 1  T ( s )

(2.11b)

Figure 2.1. A classical feedback arrangement.

Even when the selection of T(s) is based on a preliminary analysis of the closed-loop system performance with a low-order controller, there will be a considerable freedom in this selection. The selection of the required closed-loop transfer function T(s) should also satisfy the condition that the resultant controller be strictly proper or at least proper. For instance, when the plant has no open-loop zeros then T(s) should be selected such that the (poles-zeros) of T(s) must be at least greater than the plant order by one for strictly proper controller. Furthermore, there will be a cancellation of plant poles by the controller zeros as can be seen

Optimal Control Engineering with MATLAB

42

from Equation (2.11). This cancellation does not bode well for the design and for this reason such a design may be found only for naïve applications. The main problem of this cancellation is the certain degradation of the system performance as the system and/or controller parameters alter; even unstable functioning may result. Clearly, although one may use an optimized closed-loop transfer function, such a design has nothing to relate to any sort of optimized direct computation. Of course, different performances will be obtained with different selections of the optimized closed-loop transfer function. In the following example, the above procedure is carried out to obtain referential results.

Example 2.1. Consider an over damped type-zero plant which has an open-loop transfer function given by

G0 ( s) 

0.0005 (s  0.01)(s  0.05)(s  1)

It is required to design a closed-loop controlled system whose performance is defined by a minimal ITAE index for a zero steady-state step input. For strictly proper controller, a 4th-order closed-loop transfer function T(s) must be selected. From Table 2.1, we have

T (s) 

wn4 s 4  2.1wn s 3  3.4wn2 s 2  2.7 wn3 s  wn4

Hence, the controller transfer function is

GC ( s) 

2000( s  0.01)(s  0.05)(s  1) wn4 s( s3  2.1wn s 2  3.4wn2 s  2.7 wn3 )

The value of wn can be set according to the required closed-loop bandwidth, which specifies the transient response as well as the noise filtration requirements. Obviously, as the natural frequency wn increases the bandwidth increases too. The controlled system will follow the reference input faster but with less noise filtration capability. To proceed, let us set the natural frequency wn = 0.5, 1, and 2 radians/second. The system responses for unit-step input are shown in Figure 2.2. They illustrate a maximum overshoot of approximately 2% for all the three values of wn, and approximately 10.8, 5.4, and 2.8 seconds to the first overshoot respectively. The bandwidths are 0.44, 0.88, and 1.77 respectively as can be obtained from the Bode plots shown in Figure 2.3.

Controller Design Based on Parameter Optimization

43

Figure 2.2. Step responses for different values of wn (Example 2.1).

Up to this point, we shall leave the inward approach. The reader is recommended to consult the excellent treatment [6, ch.9].

Figure 2.3. Bode plots for different values of wn (Example 2.1).

Optimal Control Engineering with MATLAB

44

2.4. PARAMETER OPTIMIZATION DESIGN – OUTWARD APPROACH As it is mentioned above, another different aspect of designing an optimal control system is to assume (or to select) a specific controller structure of unknown parameters. The second step of the design procedure is to determine these parameters optimally by minimizing certain index of performance. However, to obtain satisfactory results, the selection of the controller structure and the index of performance are not straightforward. It will be always necessary to understand well the required functioning of the controlled system in relation to the plant transfer function. It is the designer‟s job to understand the customer requirements for the given plant. In this section, familiar controller structures such as PI, PD, PID, lag- network, leadnetwork, or lead-lag-network will be considered for implementing the required controller. To get a chance for computing different indices, the selection of a specific index of performance for certain given plant will not be insisted on. Nevertheless, for comparison, the reader can attempt the same plants with different indices.

Example 2.2. Consider the PID controlled system shown in Figure 2.4. Obviously, a unity- negative feedback controlled system responds very slowly to unit-step input with 0.5 steady-state error. Therefore, it is expected that a fast response will require larger gain values of the PID parameters.

Figure 2.4. A PID control system (Example 2.2).

The open-loop transfer functions is

0.0005(k 0 s  k1 s 2  k 2 ) G0 ( s )  s ( s  .01)(s  .05)(s  1) The unity-negative feedback closed-loop transfer function is

0.0005(k 0 s  k1 s 2  k 2 ) T ( s)  4 s  1.06s 3  (.0605  0.0005k1 ) s 2  0.0005(k 0  1) s  .0005k 2

Controller Design Based on Parameter Optimization

45

As a first trial of the design, the standard 4th-order polynomial of the ITAE index for zero steady-state step error systems will be used to set (optimally) the PID controller parameters. This can be done simply by equating the coefficients of equal power.

s 4  2.1 n s 3  3.4 n2 s 2  2.7 n3 s   n4  s 4  1.06s 3  (0.0605  0.0005k1 ) s 2  0.0005(k 0  1) s  0.0005k 2 Accordingly, the four unknown parameters wn, k0, k1, and k2 are calculated as

wn  1.06 / 2.1  0.5048

k1  2000(3.4wn2  0.0605)  1612 k0  2000(2.7 wn3  0.0005)  694

k2  2000wn4  130 Therefore, the transfer function of the PID becomes

130 1612(s 2  0.43s  0.0806) GC (s)  694  1612s   s s Consequently, the closed-loop transfer function is

T ( s) 

0.806(s 2  0.43s  0.0806) s 4  1.06s 3  0.867s 2  0.348s  0.065

Because of the existence of the two conjugate zeros (-0.215 ± j0.185) in the closed-loop transfer function, we do not expect a system response like that obtained in Example 2.1. However, a zero steady-state step error response is preserved. The step response is shown in Figure 2.5. The response exhibits a very high and non-acceptable overshoot of about 66%, approximately 2 and 4.3 seconds rise time and time to the first overshoot respectively, and a relatively long settling time of about 20 seconds. The task now is how to improve considerably the step response. There are many varieties; one possible technique is discussed below. It is important to note that the value of wn, which is a parameter of specifications set, is computed specifically according to the plant and the ITAE index. Therefore, there is no control to change this value for certain closed-loop requirements. This seriously limits the design. The solution to this problem is to use different controller structure, which makes the selection or the computation of wn dependent on a specific closed-loop characteristic. To damp the system further and hence to reduce the successive overshoot, let us propose a PIDD controller of the transfer function

Optimal Control Engineering with MATLAB

46

GC ( s)  K (k0  k1s 

k2 s

)( s  z )

where, the additional value of the zero z has to be found.

Figure 2.5. Step response with PID controller (Example 2.2).

Preliminary, we set K = 1. Following the above procedure, we can obtain three expressions for the controller parameters:

k 0  6800wn2  121 4200wn z  2120z k1  2000* (2.1 * wn  1.06)

2000wn4 k2  z Easily, it can be found that the zero (z) of the PIDD controller is the real root of the thirdorder algebraic equation

z 3  c2 z 2  c1 z  c0  0 where the equation constants are given in terms of wn as

c2  

3.4wn2  0.0605 2.7 wn3  0.0005 wn4 , c1  , c0   2.1wn  1.06 2.1wn  1.06 2.1wn  1.06

Controller Design Based on Parameter Optimization

47

Therefore, we can design different controllers of such a structure, for different values of wn. In our example, simulation shows that over a wide range of wn (say, 3≤ wn ≤25) the maximum percentage overshoot is reduced almost to a fixed value (approximately 46%) over the whole range, and a substantial decrease of the settling time is achieved. Figure 2.6 shows two step responses (for wn =3, and 10) , the settling time is improved significantly as wn increases. Although a reduction of 17% of the overshoot is achieved, it is still high and the controller structure has to be modified again.

Figure 2.6. Step responses for different values of wn (Example 2.2).

The next stage is to set a suitable forward gain K greater than one to reduce the overshoot and at the settling time simultaneously. For instance, we set K =5. Since the value of the ITAE index cannot be calculated analytically, it is of interest to use Simulink model to obtain a good approximation of this value. In Figure 2.7, the display reads 0.003 for the ITAE criterion. For wn = 10 radian/second, the controller is

GC ( s)  5(362800  39880s 

2515400 )( s  7.9509) s

Instead of using the gain K, an alternative usual way to improve the system response is to use a prefilter Gp(s). This prefilter eliminates the zeros of the closed-loop transfer function T(s) as shown in Figure 2.8. If the new transfer function of the closed-loop is denoted by T1(s), then it is given by the product T(s) Gp(s).

Optimal Control Engineering with MATLAB

48

T1 ( s)  T ( s)G p ( s)  Therefore, G p ( s) 

0.065 s  1.06s  0.867s 2  0.348s  0.065 4

3

12.4 , with such prefilter, the step response will be s  0.43s  0.086 2

as shown in Figure 2.9.

Figure 2.7. A Simulink model (Example 2.2).

Figure 2.8. A control system with prefilter.

The effectiveness of the prefilter on the quality of step response in terms of percentage overshoot and settling time is clear. Since the natural frequency, wn is equal to 0.505, then results almost similar to that obtained in Example 2.1 (with wn = 0.5). Finally, it is important to note that the cancellation of the controller zeros by the prefilter poles is not as serious as that done in Example 2.1 where the zeros of the controller cancel the poles of the plant. The prefilter has no contribution to the change of system parameters (sensitivity problem) and to the capability of disturbance rejection.

The control structure proposed in the previous examples is called one degree-of-freedom DOF according to, which most of the control systems are designed. However, this structure has the disadvantage that the feedback-loop properties cannot be frequency shaped independently of the reference tracking transfer function.

Controller Design Based on Parameter Optimization

49

Figure 2.9. Step response with prefilter (Example 2.2).

This means that a compromise must always be made when choosing between different time/frequency requirements. Therefore, the selection of the control and the controller structures is an important part of the design. The proposal of adding a prefilter is an example of changing the one-DOF control structure; other configurations are also possible such as with, I-PD structure. Let us now work at another example in which a PD controller of two parameters is assumed and the computation of these parameters is performed analytically.

Example 2.3. In a unity-negative feedback system, the controller GC(s) and the plant G0(s), are cascaded. The transfer functions are

G0 ( s ) 

500 , GC ( s)  K (Ts  1) s( s  16s  100)(s  5) 2

For a unit-step input, it is required to determine the parameters of the PD controller such to minimize the performance index 

J 0   y 2 (t ) dt 0

Optimal Control Engineering with MATLAB

50

By closing the loop, the system output y(t) to a unit-step input is given by the following Laplace function:

Y ( s) 

500K (Ts  1) 1 2 s  21s  180s  500( KT  1) s  500K s 4

3

According to Equations (2.4-6), the closed form of the index is given by

J0 

1 [ B0  0  B11  2b0 b1] 2K 2 

where,

 180 1 0 500K  0 500(1  KT )  21 0     det( )  0  500K 180  1   0  500(1  KT ) 21  0  K (0.82  0.695KT  0.125K 2T 2  0.11025K )  109  180 1 0 500( KT  1)  500K 500(1  KT )  21 0    0  det( )  0  500K 180  1   0 0  500(1  KT ) 21    T   K (0.2952  0.045TK )  109 K

500K  0 1  det(  0   0

500( KT  1)

0 500K  21 0  ) 0 180  1  0  500(1  KT ) 21  0.82K 2 (1  0.152439KT )  109 1

B0  25  104 K 2

, B1  25  104 K 2T 2

b0  500K

, b1  500KT

Although the index J0 has a closed form, and it will be easy to obtain the partial derivatives of J0 with respect to the parameters K and T, solving the resultant two nonlinear algebraic equations can be only carried out numerically. The MATLAB function fminsearch

Controller Design Based on Parameter Optimization

51

(which finds the minimum of a scalar function of several variables, starting at an initial estimate) can be invoked to determine the optimized values of K and T. List 2.1a, shows a MATLAB program for solving the problem in this example. The program also serves for obtaining the step responses of the control system with optimized controller and with a unitynegative feedback control only. Running the program gives the following results:

T  0.5896 , K  4.0312 , min{J 0 }  0.1698 Figure 2.10 shows the step responses. The step response of the system with the designed controller shows a faster response with a percentage overshoot of 22% and a settling time of 1.9 seconds in comparison to the over-damped response (of 2.4 seconds settling time) obtained by only closing the system with unity feedback; it might not be acceptable to have such value of overshoot.

Figure 2.10. Step responses (Example 2.3).

Example 2.4. The same problem of Example 2.3 is now solved by adopting the calculation of ISE index according to Hurwitz‟s determinants. The step response of the error and its Hurwitz‟s determinant are

E ( s) 

1 1 s 3  21s 2  180s  500  4 1  G0 ( s) s s  21s 3  180s 2  500( KT  1) s  500K

Optimal Control Engineering with MATLAB

52

0 0  21 500( KT  1) 1 180 500K 0   H a  det( ) 0 21 500( KT  1) 0    1 180 500K  0 The first row in the Hb determinant is

[bn21

 (bn22  2bn1bn3 )

(1) n2 (b12  2b0 b2 ) (1) n1 b02 ]  [ 1  81 11400  250000]

Therefore,

11400  250000 1  81 1 180  500K 0  ) H b  det( 0 21 500( KT  1)  0   180 500K  0 1 Again, the closed form of the ISE index will be some how complicated, and it is also advisable to use a numerical method to perform the optimization process. List 2.2, shows a MATLAB program which determines the two parameters K and T. According to this method, the same MATLAB function fminsearch is invoked but with a different name. The reader can check that with this program, the results are exactly the same as that obtained by program in List 2.1.

Since the MATLAB function fminsearch finds a local minimum, the selection of the initial starting point is an important issue, and many trials may be necessary before getting a suitable one. When there are many parameters, the problem becomes more troublesome. In fact, this problem is raised with all numerical optimization methods. To minimize the number of these trials, one could rely on the obtained index values which should indicate the convex property. The values of the optimized parameters themselves may also show the correctness of a specific selection. Consequently, a random number generator can be used to provide different values for initial values, and one can then select the initial vector which gives the smallest value of the index. In all cases, the initial values should give a stable response so that the index value is finite and hence the minimization process becomes meaningful. On the other hand, one may argue about the use of these closed forms to describe the index of the problem mathematically especially for highorder systems with more than one parameter. This argument is acceptable because when an efficient numerical method is used to compute the index and override some computation problems, one can get the same results. However, with the aid of MATLAB programming it is also easy to program these closed forms.

Controller Design Based on Parameter Optimization

53

Example 2.5. Consider the unity feedback system, whose open-loop transfer function is

G( s) 

40 s( s  1)(s  4)

A lag-lead compensator is designed [6] to fulfill a static velocity error constant of 10 second-1, a phase margin of 50 degrees, and a gain margin of 10 dB or more. The classical compensator design results in a lag-lead network of the transfer function

GC (s) 

(s  0.4)(s  0.2) (s  4)(s  0.02)

For our purposes, we shall assume a controller of the same structure i.e. of two zeros and two poles as

GC ( s ) 

( s  z1 )(s  z 2 ) ( s  p1 )(s  p 2 )

The Laplace transform of the error response for step input is by the rational function 4

E ( s) 

b s

i

a s

i

i 0 5

i 0

i

i

where, the denominator and nominator coefficients are given respectively as

a 0  40 p1 p 2

a1  4 p1 p2  40z1  40z 2 a2  4 p1  4 p2  5 p1 p2  40 a3  4  5 p1  5 p 2  p1 p 2

a4  5  p1  p2 a5  1 b0  4 p1 p 2

b1  4 p1  4 p2  5 p1 p2 b2  4  5 p1  5 p2  p1 p2 b3  5  p1  p 2 b4  1

Optimal Control Engineering with MATLAB

54

The ISE index will be used to determine the optimized values of the controller parameters. List 2.3a, shows a MATLAB program in which the method of Hurwitz‟s determinants is programmed, and the MATLAB function fminsearch is also invoked to perform the minimization process. The initial values are selected so that the initial error step response is stable (the reader can check that if the initial vector is [0.5 0.5 0.5 0.5], then an unstable response is observed, and the minimization process is no longer correct). In addition, in the program the result obtained in [6] is included for making a comparison. The designed controller transfer function is

GC ( s ) 

( s  0.34535)(s  0.34538) ( s  0.3454) 2  ( s  3.8168)(s  0.17464) ( s  3.8168)(s  0.17464)

The step responses of our closed-loop system, and the one obtained in the reference are shown in Figure 2.11. As it can be seen, the optimized system has very fast settling time with almost the same rise time and overshoots. As it is expected this will reduce the relative stability of the system. The MATLAB function allmargin for both systems indicates a loss of gain margin and phase margin of about the 2 dB and 36 degree respectively.

Figure 2.11. Step responses of optimized and classical designs (Example 2.5).

The same results are obtained if the error values of the unit-step response are squared and numerically integrated over the time interval 0  t  20 seconds. List 2.3b, shows the m. program (Four_PC(v)), which is used instead of (Four_P(v)) in List 2.3a to compute the square error numerically by Simpson‟s method. The fixed step of integration and the odd number of data which are required by the used numerical method could be obtained by the MATLAB function step. The conclusion here, by the parameter optimization technique we can synthesize a controller of a given structure directly without any sort of trial-and-error procedures. The

Controller Design Based on Parameter Optimization

55

optimized controlled system can have similar results as that designed by, for example, usual root locus, Bode plots, or other frequency techniques The reader may wish to know the system response for different indices using such numerical computation. This can be simply performed by changing the statement of computing the variable integd in the program of List 2.3b.

2.5. LIMITATIONS OF PARAMETER OPTIMIZATION For designing a proper stable controller of as simple as possible structure, there are several limitations in applying the parameter optimization approach. Even when we exclude the problems concerning the numerical approach of finding the minimum of the index such as a proper guessing of the initial values or the computation accuracy, there will be still an important question about the convergence of the solution to realistic values. The convexity of the closed-form function in a convex region R may not indicate that the minimization process will converge to the required answer. Mathematically, the minimization process may converge to positive real or negative real values (the solution of nonlinear equations obtained from the minimization may even converge to complex values), while the controller parameters have to be positive real values (for strictly non-minimum phase control). To perceive such serious limitation, the next example is introduced.

Example 2.6. Consider a plant G(s) and a suggested cascade GC(s) controller in a unity-negative feedback system

G(s) 

1 s (T p s  1)

,

GC ( s ) 

K , T p , Tc , K  0 Tc s  1

It is required to determine the controller parameters Tc and K such that ISE index is minimal. The Laplace transform of the error step response is

E (s) 

Tc T p s 2  (Tc  T p ) s  1 Tc T p s 3  (Tc  T p ) s 2  s  K

Using Hurwitz‟s determinants, the closed form of the index can be found in terms of system and controller parameters as

Optimal Control Engineering with MATLAB

56

(Tc  T p ) 2 1 J0   2 K 2(Tc  T p  KTc T p ) Then, the two partial derivatives of J0 with respect to the two parameters TC and K are determined and are set equal to zero.

J 0 (Tc  T p )[(Tc  T p )( KT p  1)  2 KTc T p ]  0 Tc 2(Tc  T p  KTc T p ) 2 Tc T p (Tc  T p ) 2 J 0 1   0 K 2 K 2 2(Tc  T p  KTc T p ) 2 Solving these two equations for Tc and K we get the following two coupled equations:

Tc  K

T p ( KT p  1) KT p  1 Tc  T p Tc T p  Tc T p (Tc  T p )

With simple algebraic manipulation, we can get one equation for Tc

Tc  T p Tc  T p  0 3

3

Since the plant time constant Tp is a positive real value, then the above equation has no real positive root and hence the value of the controller time constant Tc is not realizable. In other words, such parameter optimization cannot be performed. To proceed let us assume that Tp = 0.04, and we choose an arbitrary value of Tc say 0.05. Now, we want only to optimize the controller gain K. Then, the above solution gives

K

0.05  0.04 (0.04)(0.04)  (0.04)(0.05) (0.05  0.04)

 14.938

Substituting these values {0.04, 0.05, and 14.938} in the partial derivatives, we get

J 0 0 , K

J 0  0.987, not both equal to zero Tc

Controller Design Based on Parameter Optimization

57

Simple algebraic analysis can show that as Tc increases both of the gain K and the partial derivative with respect to Tc decrease monotonically. Therefore, the choice of a specific value of Tc becomes a performance issue. Alternatively, for given Tp, one may try to determine the values of K and Tc, which give the smallest possible real value of the function defined by

J m  (Tc 

T p ( KT p  1) KT p  1

) 2  (K 

Tc  T p Tc T p  Tc T p (Tc  T p )

)2

For Tp = 0.04, the smallest positive value of Jm equal to 0.0404, when the controller parameters are K = 0.7792, and Tc = 0.0524. With these parameters, the step response of the closed-loop system is shown in Figure 2.12. As it is known for control engineers, the selection of specific controller structure to fulfill certain closed-loop specifications can be made mostly according to the known facts in control theory. However, as one requires determining the controller parameters by minimizing certain index of performance the problem of selecting a specific index is encountered. From the above examples, we notice that the minimization of ISE gives relatively high overshoot and fast response for step input i.e. less relative stable performance. We may want to have a similar value of overshoot (without using prefiltering technique as used in Example 2.2) via minimizing another index of performance.

Figure 2.12. Step response (Example 2.6).

Generally, there are no longer confined guidelines of doing so. Let us now work again the previous example but with different index of performance.

Optimal Control Engineering with MATLAB

58

Example 2.7. It is required to determine the controller parameters Tc and K of the control system in Example 2.6 such that ITE index is minimal. A closed form of ITE index is given by 

J   te(t )dt  lim(1) s 0

0

E (s) 

d d E ( s)  lim(1) E ( s) s   ds ds

Tc T p s 2  (Tc  T p ) s  1 Tc T p s 3  (Tc  T p ) s 2  s  K

Therefore the closed form of the ITE index is

J

1 (Tc  T p )  K K2

Then we have

(Tc  T p ) J 2 2  3  0  K  2 K Tc  T p K K Note that there is no meaning of optimizing with respect to Tc, and hence it is similar to what we obtained in Example 2.6. However, in this case, the value of Tc contributes clearly in controlling the step response. Although this one degree-of-freedom case (shown in the previous example) gives design flexibility, the solution is no longer optimal. Moreover, with a three parameters controllers and when the minimization of ITE index has more than one degree of freedom, then even this flexibility will disappear. The convergence to unreliable values includes not only the negative values, but also the very large positive values, which could cause actuator saturation. Therefore, searching about a specific index, which fulfills a certain set of requirements and overpasses all convergence problem, is a challenging issue. The next example elucidates such limitations.

Example 2.8. Consider the following plant and a PID controller in unity-negative feedback system

G( s) 

1 2 s  3.6s  9

,

GC (s) 

k1 s 2  k 2 s  k 3 s

Controller Design Based on Parameter Optimization

59

It is required to determine the PID parameters k1, k2, and k3 such that ISE index is minimal. A closed form of the ISE index is

J

k 3 k 2  3.96k 3  81k1  291.6 2(k1 k 2 k 3  9k1 k 2  3.6k 3 k 2  32.4k 3  k 32 )

As it can be deduced from the above expression, the value of the ISE index is decreasing toward its minimal value as the three gains increase. Therefore, the minimization process continues while the three parameters increase without bounds. However, if a small tolerance value is requested to end the minimization process, then we will obtain very large values of the three parameters. For instance, in this example, if a numerical integration method is used to compute the index to a finite sufficiently large time (say 5 seconds), then the gain values are of order 109, and the minimal value of the index is of order 10-3. Obviously, such large values will cause actuator saturation. Of course, if one postulates a large value for the index, then the optimization meaning will be lost. If the ITSE index is proposed, similar ill-conditioning results will be obtained. Let us propose the ITE index as we did in Example 2.7. We can show that the following expression is valid for computing the index

J 

3.6 9k 2  81  k3 k 32

Hence the index derivative can be performed only for k3,

J 3.6 2(9k 2  81) 0 2   0  k 3  5k 2  45 k 3 k3 k 33 As we said before, there will be two degree-of-freedoms in this case. This is the choice of the two parameters k2 and k1. Hence, we contrive to do some trials to achieve a required response. Let us now propose the ITAEER index to determine the PID parameters. The numerical calculation of the error derivative is performed with 5 points of Newton‟s polynomial. Then the total expression of the integrand is calculated numerically using Simpson‟s method (one can invoke the MATLAB function Trapz which gives slightly different results). To avoid illconditioning, a check for the existence of a positive pole in the characteristic equation is carried out. If such an occurrence is detected, then this will indicate the incorrect direction of searching a minimum. List 2.4 shows the corresponding MATLAB program. Without any detection of ill-conditioning cases, the running of the program gives the following results:

k1  151.9 , k 2  546.7 , k 3  1367 , J  0.0089

60

Optimal Control Engineering with MATLAB

Let us assume that these values are still large and saturation may occur. A simple solution is to divide these values by a factor q. The step responses for different values of q are shown in Figure 2.13. When q = 15, we get the largest settling time which is equal to 0.39 seconds. With q > 1 all such solutions are suboptimal. Therefore, the designer may accept the PID controller whose transfer function is

U ( s ) 10.127s 2  36.447s  91.133  E (s) s

Figure 2.13. Step responses for different values of factor q (Example 2.8).

For the same considered plant, a PID controller is designed using the root locus technique [1]. The control purpose is that the control system is capable of reduce sufficiently a step input disturbance. The control system with the input disturbance is as shown in Figure 2.14.

Figure 2.14. A PID control system.

It is of interest to compare results. Figures 2.15a and b show respectively the unit-step responses of both controllers, and the system responses to unit-step disturbance input.

Controller Design Based on Parameter Optimization

61

The optimized design has a dead-beat response to a unit-step reference input with approximately 0.4 seconds settling time. This is compared to 7% overshoot and 0.8 seconds of the root locus design. For unit-step disturbance input, the optimal response settles faster, and has less magnitude. Therefore, the overall performance of the optimized design is more acceptable than the root locus design. In all previous examples, the controller is connected in the forward path. However, the controller may be connected in the feedback path or even two different controller structures may be suggested to be connected one per a path. In classical design methods such as the root locus or frequency loop shaping of Bode plots or Nyquist plot … etc., the selection of one of these configurations may be performed depending upon a set of thumb of rules. With parameter optimization methods, such a thumb of rules will not be so useful because of the causality missing between the proper selection of the method of optimization and the result of the optimization process.

Figure 2.15. Step responses with reference input (above) and disturbance (below) (Example 2.8).

For example, during the practical implementation of a PID controller, we may wish to operate the derivative action only in the feedback path so that differentiation occurs only on the feedback signal but not on the reference signal. The control scheme arranged in this way is called PI-D. For other reasons, we may also wish to put the proportional term in the feedback path to make what is called, I-PD scheme. Another example, of changing the control system configuration is found when we wish to reduce the effect of disturbance and/or noise. Although in all such cases, the characteristic equation is kept unchanged the signals in feed forward or feedback paths will be altered. Consequently, the output or the error responses are varied from one configuration to another.

62

Optimal Control Engineering with MATLAB

In the next example, we shall consider an I-PD scheme, which is proposed to control the plant in Example 2.8. Based on Example 2.8 results, the idea is to discuss some difficulties or perhaps ambiguities in applying parameter optimization.

Figure 2.16. I-PD control system (Example 2.9).

Example 2.9. Consider the I-PD control system shown in Figure 2.16. In this control system, the output and error unit-step responses are given by the following Laplace transforms:

Y (s) 

k3 1 2 s  (k1  3.6) s  (k 2  9) s  k 3 s

E ( s) 

3

s 2  (k1  3.6) s  (k 2  9) s 3  (k1  3.6) s 2  (k 2  9) s  k 3

When the adopted values of k‟s obtained in Example 2.8 are used, then the output and error unit-step responses together with that obtained in Example 2.8 are as shown in Figure 2.17. As it can be seen both the output and error responses are changed. For instance, the maximum overshoot is increased by 9% and the settling time is also increased to about 2.3 seconds. To adjust again such degradation, we may think about re-optimizing these values according to the new error function by minimizing the same index ITAEER. For consistency, we shall adopt the same initial conditions for starting the minimization process and the same numerical method with the same computation time interval, i.e. to use the same program shown in List 2.3. Ill-conditioning of minimization, which is defined by a positive pole in the characteristic equation, is detected during the minimization process and hence searching a minimal value is stopped. Since mathematically the search for a minimum can be continued in spite of this illconditioning then we hope that a convergence to realistic values might be obtained if the

Controller Design Based on Parameter Optimization

63

condition of positive pole occurrence is discarded. Doing so, the searching for a minimum continues, and it converges to the following results:

k1  8.3003 , k 2  105.6603 , k 3  331.2532 , J  0.4304

Figure 2.17. Output and error step responses (Example 2.8 and Example 2.9).

The minimal value of the index is very large in comparison to that obtained for the PID controller (see Example 2.8). The step responses of both I-PD and PID controllers are shown in Figure 2.18. The control system with I-PD controller becomes more sluggish. The settling time becomes greater than one second. Therefore, we can say that for this particular case, the parameter optimization fails to improve the system response.

One may argue about changing the initial starting values or changing the computation time or even the index, or to use more than one degree-of-freedom controller. For instance, the reader may examine the use of two PID controllers in the inner and outer loops. However, making such trials bereave the procedural privilege of this approach. We can finally conclude that the parameter optimization approach depends on the considered system, and the control engineer has to explore several alternatives. The exact solutions can be obtained only for some cases. Furthermore, for nonlinear control systems, it is not so attractive approach as with linear control systems. In spite of these drawbacks, this design approach is widely used in the field due to the availability of various methods of optimization. It is compatible with the linear control theory. Moreover, this design approach keeps pace with the new methods of optimization.

64

Optimal Control Engineering with MATLAB

Figure 2.18. Step responses of optimized I_PD and PID controllers (Example 2.9).

2.6. PARAMETER OPTIMIZATION WITH CONSTRAINTS As it is noticed from the previous section, there are many cases where the minimization process cannot specify one or more of the controller parameters. Moreover, in other cases, these parameters may have large values, which can easily saturate the system actuators. This includes the zero values of some parameters and/or unbounded values of the other. In any case, the assignment of practical values changes the optimal solution to a suboptimal one; as it was done in previous examples (viz. Example 2.7 and Example 2.8). Such problems can be solved either by changing the index and/or specifying certain physically meaningful constraints to be satisfied. On the other hand, in some applications, the constraints are inherent parts of the control system. For example, in electromechanical systems, the electric current or the shaft velocity should not exceed certain values. In other examples, mechanical constraints enforce the angular or the linear displacement motion to be stopped at a certain position to prevent damages. The bandwidth of the closed-loop system, the H2 or H∞ norms of the closed-loop, a signal at some place in the control system, a required energy defined by the integral square of the input control u(t), etc., are often bounded to specific values. In this section, we would like to highlight the solution of the parameter optimization problem when there is a specific constraint on the input energy. Such constraints are often used to prevent the saturation of actuators. For constant value M, it has the form 

J u   u 2 (t )dt  M 0

(2.12a)

Controller Design Based on Parameter Optimization

65

Therefore, the parameter optimization should be performed to minimize a given index and at the meantime satisfying the energy constraints for the steady-state case. Mathematically, it is known that the minimization process can be still carried out by connecting the index and the constraint multiplied by Lagrange's multiplier in one single new index, i.e. 

J n  J 0    u 2 (t )dt

(2.12b)

0

where J0 can be any of the error indices shown previously, which should be given in a closed form in terms of the system and controller parameters. The optimization is carried out by solving simultaneously the equations system given by

J n J 0 J    u  0 , i  1,2..., P pi pi pi

(2.13)

J u ( pi )  M where pi are the ith controller parameters, and P is number of parameters. Therefore, rather than determining some control u(t), from a class of control functions, which satisfy the constraint, to minimize Jo, an equivalent problem is solved. It is to choose u(t) from a class of control functions, which minimize Jo for some value of λ, and satisfy the constraint of Equation (2.11). Note that the value of λ will be found from Equations (2.13).

Example 2.10. Consider the unity-negative feedback control system shown in Figure 2.19.

Figure 2.19. A control system (Example 2.10).

The open-loop transfer function is

G0 (s) 

K s(0.5s  1)

It is required to specify the parameter K such that the ISE index is minimal.

Optimal Control Engineering with MATLAB

66

Similarly, as we did before, the index can be determined to have the following closedform,

J

1 1  4 2K

J 1 0  K  K 2K 2

Therefore, the minimization of such an index has no practical meaning. Instead we shall minimize the ISE index and at the same time satisfying the following constraint 

J u   u 2 (t ) dt  10 0



J n   [e 2 (t )  u 2 (t )]dt 0

Using the Hurwitz‟s determinants, the closed forms are obtained.

1 1 J0   4 2K K J u   10 2 1 1 K Jn    4 2K 2 Although in this simple example K = 20 is found easily, we proceed to determine the partial derivative of Jn. J n 1    0 2 K 2 2K 1   2  0.0025 K The minimal value of Jn is 0.03.

Example 2.11. Consider the control system shown in Figure 2.20. It is required to determine the controller parameters K and T such as to minimize the ISE index and at the same time satisfying the condition

Controller Design Based on Parameter Optimization

67



J v   v 2 (t ) dt  25 0

The actuating error e(t) and the velocity v(t) signals have the following Laplace transform to unit-step reference input r(t)

ss((Ts Ts  1 1)) E ( s )  3 2 E ( s)  Ts 3  s 2  5KTs  K Ts  s  5KTs  K K K ((5 5Ts Ts  1 1)) V V (( ss))   Ts 33  s 22  5KTs  K Ts  s  5KTs  K

Figure 2.20. A control system (Example 2.11).

The closed forms of the ISE index and the constraint are given respectively by

5KT 22  1 JJ 0   5KT  1 0 8KT KT 8 25KT 22  1 JJ v   25KT  1   25 25 v 8T T 8 Therefore a combined performance index will be

J n  J 0  J v 

5KT 2  1  25K 2T 2  K 8KT

The partial derivatives of Jn with respect to the controller parameters K and T are

J n 1 KT (5T 2  50KT 2   ) (5KT 2  1  25K 2T 2  K )T  {  } 0 K 8 K 2T 2 K 2T 2 J n 1 KT (10KT  50K 2T ) (5KT 2  1  25K 2T 2  K ) K  {  } 0 T 8 K 2T 2 K 2T 2

68

Optimal Control Engineering with MATLAB The two derivatives give respectively

25K 2T 2  1 25K 2T 2  5KT 2  K  1 Solving for K and T in terms of λ, we get

K T

1 5

 5

By substituting these values in the given constraint, the value of the optimum λ can be then calculated

25

1 5 8

(

 5



)2  1  25    0.00131

5

Finally, the optimum controller parameters and the minimal ISE values are obtained.

K  341.64, T  0.0162, J 0  0.033 Figure 2.21a shows the step output response, which has a relatively high overshoot, 36%, and it is settled within 0.4 seconds. Figure 2.21b shows the signal v(t). The value of the constraint Jv can be computed numerically with sufficiently small step of integration to indicate the validity of the solution. Doing so we will obtain a value of 25, i.e. the given constraint is satisfied over the whole time of control. The reader can make a Simulink model to observe that for any simulation time, the value of the Jv will not exceed 25. Furthermore, it can be shown that as the constraint Jv becomes tighter (e.g. Jv ≤ 5), then the response will be more sluggish (settled within two seconds) and as a result of this the index Jo will be approximately five times greater. The same limitations discussed above existing in this type of optimization. Moreover, the analytical solution of the (P+1) nonlinear algebraic equations is not always easy to obtain especially when the number of the unknown parameters exceeds two. Even when numerical methods are adopted, the selection of the initial start values is the most difficult problem, and perhaps one may be lucky if he or she could find the solution in a finite number of trials. Furthermore, even if the numerical solution is searched directly through Equation (2.14) then the same numerical problems arise. The MATLAB optimization toolbox offers the commands to solve the problem by utilizing the MATLAB function, fmincon (The minimization with constraints function). This

Controller Design Based on Parameter Optimization

69

command search for a local minimum of multi-objective function subjected to equality and/or inequality constrains; many helpful facilities are also included [7]. For the problem considered in Example 2.11, the program in List 2.5 is written to determine the two parameters K, and T directly by specifying Jo as the function to be minimized and Jv as a constrain. Running the above program gives after 24 iterations (and 100 times function counting) the same results that are obtained analytically. It is worth mentioning that the initial vector x0 is not arbitrary. It should be properly selected in order to terminate accurately at a local minimum. For example, while x0 = [1.5 0.5] terminate to slightly different results, the initial vectors, x0 = [0 0] or x0 = [2 1] are both incorrect selections.

(a)

(b) Figure 2.21. a) Output response y(t) b) The signal v(t) (Example 2.11).

70

Optimal Control Engineering with MATLAB

Therefore, although, the parameter optimization approach seems to be valuable and easy to carry out, the convergence problem could be serious in certain circumstances. Moreover, practically, it leads often to suboptimal solution. More effective way will be described in chapter five.

2.7. CONTROL VECTOR PARAMETERIZATION In parameter optimization approach, we worry about the controller structure and not about the output of the controller. Therefore, an alternative approach is to carry out the parameter optimization by the parameterization of the control variable u(t) as a function of time. This approach was first announced and used in [8, 9]. It was applied to solve chemical engineering problems. Since mathematically the function u(t) can be approximated by a finite series of orthogonal time functions then we can write i n

u (t )   ai f i (t )

(2.14)

i 0

where ai are unknown coefficients and fi(t) are orthogonal functions such as Laguerre functions, Hermite polynomials, Chebyshev polynomials, etc. The coefficients (parameters) are then found for minimal value of the proposed index. The choice of the number of functions should be based on a prior knowledge of the u(t) function itself in order to have a close approximation to the optimum with only a small number of terms. In fact, such a condition is perhaps the major drawback of the method. Although, Chebyshev polynomials require equal spaced data, they will be adopted here for approximating the optimal control input. The Chebyshev polynomials of a first kind have the forms [10]

T0 ( x)  1 T1 ( x)  x T2 ( x)  2 x 2  1

(2.15)

: Tn ( x)  2 xTn 1 ( x)  Tn  2 ( x), n  2,3,... where the x lies in the range [-1, +1]. To convert from a different range of the independent variable t, say, in the range [a, b], then one can use the following conversation formula:

x

2t  (a  b) (b  a)

(2.16)

Controller Design Based on Parameter Optimization

71

For instance, if our problem is to perform the optimization between t0 = 0 and tf = 3 seconds then the used Chebyshev polynomials will have the forms

T0 (t )  1 2 t 1 3 2 T2 (t )  2( t  1) 2  1 ... etc 3 T1 (t ) 

(2.17)

In our problem, this conversion is essential mainly because the initial time is often zero. The next example illustrates the parameterization idea of control variable. This example is of the type of the optimal control problems that will be considered later in this book.

Example 2.12. Consider the nonlinear first-order system described by the state equation

x   2 x 2  u , x(0)  1 It is required to determine the control variable u(t), which transfers the system to an arbitrary final state value while minimizing the index 2

J (u )   (0.75 x 2  0.5 u 2 ) dt 0

We will assume that the control variable has the form

u(t )  a0  a1t  a2 t 2 The state equation will be solved numerically by the MATLAB function ode45. The minimization of the index will be again performed by the MATLAB function fminsearch. The index value is computed numerically using Simpson‟s method. The List 2.6 shows a complete MATLAB program for solving the problem. The results are:

u (t )  0.5711 0.6422t  0.1965t 2 The final state value x(2) is 0.094 and the minimal value of the index is 0.1376. Figure 2.22 shows, the optimal trajectories over the time interval [0, 2].

Optimal Control Engineering with MATLAB

72

Figure 2.22. Optimal state trajectories and control (Example 2.12).

Since the system state can have any real value at the final time, and the minimal value of the index is grater than zero. Therefore, one may inquire about the existence of another control variable function, which gives a minimal value smaller than the obtained value. Simulation shows that increasing the number of terms in the control variable function reduces the minimal value; for instance, with 4, 5 terms, the minimal values are 0.1357 and 0.1336 respectively. These figures indicate that the chosen 3 terms control variable function is sufficient. When the problem is constrained by some terminal state conditions, the minimization process should be constrained also by these conditions, otherwise this method will fail to solve the problem. This can be simply performed by adding these terminal constraints with suitable (adjustable) weights to the index. The next example illustrates the use of Chebyshev polynomials to approximate the variable control when the terminal state constraints are also required.

Example 2.13. For the double integral control system described by

x1  x2 x2  u

, x1 (0)  0 , x2 (0)  0

We want to determine the control variable u(t) such that in the time interval [0, 2], the system has to reach the terminal values x1(2) = 1, and x2(2) = 0, and minimizing the index

Controller Design Based on Parameter Optimization

J (u ) 

73

2

1 ( x12  u 2 ) dt  20

Let us assume that an approximated control variable u(t) is given by

u (t )  a 0T0 (t  1)  a1T1 (t  1)  a 2T2 (t  1) where the argument (t-1) renders Chebyshev polynomials orthogonal in the interval of the example [0, 2]. Therefore, the control variable and the system states are

u (t )  (a 0  a1  a 2 )  (a1  4a 2 )t  2a 2 t 2 1 2 x 2 (t )  (a 0  a1  a 2 )t  (a1  4a 2 )t 2  a 2 t 3 2 3 1 1 1 x1 (t )  (a 0  a1  a 2 )t 2  (a1  4a 2 )t 3  a 2 t 4 2 6 6 The index should be modified to take into account the terminal conditions. The new index form will be

J  [k1 ( x1 (t )  x1 (2)) 2  k 2 ( x2 (t )  x2 (2)) 2 ] |t t

2

f

0.5 ( x12  u 2 )dt 0

where k1 and k2 are two penalty weighting coefficients. During the iterations of minimizing the modified index, the states will have different values than the terminal values. The differences are weighted equally by k1 = k2= 1000 to minimize these differences at the terminal time. These weighting values give the following results:

u(t )  0.1367T0  1.4976T1  0.411T2  0.951 0.147t  0.823t 2 The index minimal value is 0.9182, and the terminal states are reached. If closer results to the required terminal states are required, then the values of k1 and k2 must be increased further (perfect terminal matching requires that k1→ ∞, and k2→ ∞). Figure 2.23 shows the optimized trajectories of the control variable and the system states. We would like to exploit this example to illustrate how to implement the optimal time control as usually found in control systems as a transfer function. The optimal open-loop linear system is obtained by designing a controller which would convert a unit-step input to u*(t), obtained above. The resulting system is shown in Figure 2.24a, where,

G0 (s) 

1 0.951s 2  0.147s  1.646 1 , G ( s )  , R( s )  C 2 2 s s s

74

Optimal Control Engineering with MATLAB

Figure 2.23. Optimal control and state trajectories (Example 2.13).

The transfer function of the optimal controller in a unity-negative feedback closed-loop system, which is shown in Figure 2.24b, is obtained easily from

Gc f ( s ) 

GC 0.951s 4  0.147s 3  1.646s 2  4 1  GC G0 s  0.951s 2  0.147s  1.646

(a)

(b) Figure 2.24. Optimal control systems a) open-loop, b) closed-loop (Example 2.13).

Controller Design Based on Parameter Optimization

75

Of course, both systems have the same U*(s) and the same unit-step responses y(t), where y(2) = 1 as required. As a final note, let us explore the case where it is also required to satisfy an energy constraint in the form 2

0 e

2

(t ) dt  M

where e(t) = r(t) – y(t). In this example, the value of this integral is 0.82, so the constraint is satisfied for M = 1, Therefore, no additional efforts are required to satisfy this constraint. However, let this constraint is replaced by 2

0 u

2

(t ) dt  M

Then its value will be 1.615 which is greater than M = 1, and in this case the designer has to find an adequate other method to solve the problem; the technique adopted in Example 2.11 is not the only way.

2.8. PARAMETERS OPTIMIZATION VIA GENETIC ALGORITHMS In this last section, it is worthwhile to include a different optimization approach based on one of the evolutionary techniques, namely the genetic algorithm GA. Genetic algorithms were invented and developed by John Holland and his colleagues in the 1970s. Later, the excellent book of Goldberg [11] brought the subject to the application area. Porter and Jones [12] proposed the GA as an alternative means of auto tuning of PID controller. Such use of GA became attractive with different controller structures [13]. Genetic algorithms are known to be robust and have been enjoying increasing popularity in the field of numerical optimization. These search algorithms are inspired by natural selection and natural genetic (emulate natural genetic operations), and it is more probable to converge to a global solution. One of GA key features is that it searches from population of points rather than from a single point. Distinct point of interest in using GA for optimizing controller parameters is that it uses objective function information in more flexible and general forms instead of derivatives. Moreover, it is equally applied to time-invariant, timevarying linear or nonlinear systems. However, besides the known fact of the large computation time, GA has two additional drawbacks, which are well understood as they affect the optimization process. Some objective functions are GA-deception functions, where bad chromosomes can be generated by combining good building blocks, which cause the failure of the building block hypothesis. A small population size may cause what is called the genetic drift phenomenon, and therefore, obtaining a global point may not be guaranteed. It is not the task of this book to discuss GA‟s operations and features. Instead, we utilize the GA MATLAB toolbox to optimize the parameters of fixed structure controllers such as to

Optimal Control Engineering with MATLAB

76

satisfy given performance index. However, in short, we say that GA has three main operations. They are the Selection, Crossover, and Mutation operations. A simple genetic algorithm iterates as follows: 1. 2. 3. 4. 5.

6. 7. 8.

Start with randomly generator population of n chromosomes. Calculate the fitness of each chromosome in the population. Repeat the following steps until n offspring have been born. Select a pair of parent chromosomes from the current population based on the increasing of the fitness function. With crossover probability, crossover the pair at a randomly chosen point to form two offspring; if no crossover occurs the two offspring are exact copies of their respective parents. Mutate the two offspring at each locus with certain small probability and place the resulting chromosome in new population. Replace the current population with new one. Go to step 2.

Readers who are familiar with GA‟s theory can simply make use of the MATLAB toolbox guide to design their controllers, otherwise it is recommended first to refer to the special literatures, for example, [15]. In particular, two MATLAB functions are sufficient to perform the task. The MATLAB function ga finds the minimum of the required performance index. The genetic algorithm parameters can be set the MATLAB function gaoptimset. To illustrate, tuning a PID controller via GA technique is introduced in the next example. A comprehensive treatment of PID tuning based GA‟s including different PID structures can be found in [15].

Example 2.14. In a unity-negative control system, the plant and the controller are described by the transfer functions

G0 ( s ) 

5

( s  3)(s  3s  3) K GC ( s )  K p  i  K d s s 2

First let us tune the PID controller such as to minimize the ITAE error. The problem is solved using the MATLAB genetic toolbox as shown in List 2.7. The GA utilizes a default set of parameters except the number of generations and the population size (300 and 30 respectively). The m. file named genfun2 computes the index ITAE to be minimized by collecting the data of the step error response in integration step of 0.01 in the interval [0, 20], and performing the usual procedure of numerical integration. The reader will notice that the execution time of the MATLAB program is somehow long as compared to previous programs. In fact, this is the main drawback of utilizing GA in real time control.

Controller Design Based on Parameter Optimization

77

The results are as follows:   

The tuning values of PID parameters are:

K

p

Ki



K d  6.6858 5.3699 2.9356

The minimum value of the index is 0.1813. The transfer function of the closed-loop controlled system is

14.68s 2  33.43s  26.85 Gcl (s)  4 s  6s 3  26.68s 2  42.43s  26.85 

The step response is shown in Figure (2.25), where an overshoot of 9% and settling time of 2.13 seconds are noticed.

Next let us require minimum ITAE, and moreover, minimum overshoot. It is easy to modify the m. file genfun2 to perform this requirement; the overshoot can be computed by the MATLAB functions abs and min (overshoot = abs (min (e))), and the total cost function becomes 20

J m  w1 (overshoot)  w2  t | e | dt 0

where w1 and w2 are appreciated weights. For instance, for w1 = w2 = 1, we get

K

p

Ki



K d  5.7539 4.5511 2.7848

J min  0.2376 overshoot  4% setteling  2.24 sec. Gcl ( s ) 

13.92s 2  28.77s  22.76 s 4  6 s 3  25.92s 2  37.77s  22.76

The step response is shown in Figure 2.26. Different values of the weights w1 and w2 result in a competition between the overshoot and the settling time, for example, w1 = 0.8 and w2 = 0.2 result in an overshoot of less than 1%, and 2.29 seconds, settling time. Finally, using the closed-loop 4th-order characteristic equation listed in Table 2.1 one can determine the PID parameters (see Example 2.2) which give 39% overshoot and a settling time of 2.96 seconds. Therefore, the powerful of GA as a minimization approach is so clear.

78

Optimal Control Engineering with MATLAB

Figure 2.25. Step response for minimum ITAE (Example 2.14).

Figure 2.26. Step response for minimal Jm (Example 2.14).

It is worth mentioning some important points. First the MATLAB function gaoptimset should define the adequate set of GA parameters. In special problems such as nonlinear systems, it may be required to set only a specific set of parameters. Second, although in some cases it may not be required to set the range values of the controller parameters, in other cases such a setting is essential. Furthermore, the execution time differs from one problem to another; actually, it depends not only on the number of generations but also on the other GA parameters and on the problem itself. Finally, due to the random nature of the GA the results may vary slightly each time the genetic algorithm is executed. The MATLAB offers the

Controller Design Based on Parameter Optimization

79

possibility of reproducing the same results by restarting the rand and randn generators. Nevertheless, the randomness nature could be useful for avoiding local minimum traps.

PROBLEMS 2.1. Consider a ship roll stabilization system shown in Figure 2.27.

Figure 2.27. Ship roll stabilization system (Prob. 2.1).

The stabilizing fin actuator and the ship roll dynamic are respectively

Ga ( s ) 

1 1 , Gs ( s )  2 s 1 s  0.7 s  2

If root locus technique is used to design the stabilizing controller, it is found [16] that a PIDD controller of transfer function

k (s  z )(s 2  as  b) GC  , k  10.2, z  2, a  4, b  8 s achieves less than 25% overshoot and less than 2 seconds settling time. i. Determine the k, z, a, b parameters such that the ITAE performance index is minimal, where e(t) = ζd(t) – ζa(t). ii. Improve the transient and steady-state characteristics by a suitable forward gain K. iii. Evaluate the optimum value of ITAE index and compare with that obtained with the controller designed in the reference. 2.2. In a unity-negative feedback system, the controller GC(s), and the plant G0(s) are cascaded. The transfer functions are

G0 ( s) 

10 sa , GC ( s)  s( s  1)(s  5) s  a

80

Optimal Control Engineering with MATLAB 

0

The performance index is J 0  e (t )dt to be minimal for unit-step input. 2

i. Determine the controller parameters and state, whether it is a phase-lag or a phase-lead network. ii. Set up a Simulink model for the controlled system, and determine the maximum overshoot and the settling time. 2.3. Consider the control system shown in Figure 2.28, and a performance index of the form 

J 1   [qe 2 (t )  u 2 (t )] dt 0

If r(t) is a unit-step reference input and q = 25, show that

J1 

12.5( K 3  465K 2  168K  1182) ( K 2  465K  139)

In addition, for minimal J1, the optimum value of the feedback gain K is 1.29. Simulate the system to determine the unit-step response. Give comments concerning the effect of q value on the speed of response.

Figure 2.28. A Control system (Prob. 2.3).

2.4. Consider the block diagram shown in Figure 2.29 of an electrically propelled car. The speed of the car is to be controlled by actuating two AC motors of a transfer function

G( s) 

6 ( s  0.5) 2

Controller Design Based on Parameter Optimization

81

Figure 2.29. A car control system (Prob. 2.4).

The tacho-generator transfer function is assumed to be equal to 0.85. Design a PID controller based on the ITAE performance criterion such that the peak overshoot is less than 2%. Hint: if it is required, use a preamplifier. Plot the transient response for a step command. 2.5. A servomechanism has an open-loop transfer function

G0 (s) 

200 s(s  2)(s  10)

The proposed control system is a unity-negative feedback with feed forward network of a transfer function

GC ( s ) 

1  T1 s 1  T2 s

i. Show that the performance index ITE can not be used for determining the controller parameters. ii. Based on the minimization of the ITSE, determine the time constants T1 and T2. Simulate to obtain a step response. iii. Suggest an additional control to improve the system performance such that a PM > 45○, and GM > 20db are achieved. 2.6. Resolve the problem in Example 2.8 but with ISEER index. 2.7. Consider the control system shown in Figure 2.30.

Figure 2.30. A Control system (Prob. 2.7).

82

Optimal Control Engineering with MATLAB

i. Show that neither the ISE nor the ITE indexes can be used to determine the controller parameters, while the index ITSER dose. ii. For the controlled system designed in i) above, If a 0.03 forward gain is introduced, determine whether the system satisfies the constraint, 

S   u 2 (t ) dt  50 0

Use Simulink model to state your answer. 2.8. A simplified yaw damper mode of a general aviation aircraft is represented by the control system shown in Figure (2.31). Determine the controller gain K which will result in a minimum ISE and fulfilling the constraint δR ≤ 10o for a unit-step command rc. Hint: The E(s) transfer function of Equation (2.9) and the ISE form should be modified to

E ( s) 



A 2  Em ( s), J 0   em (t ) dt s 0

In addition, the maximum output value of the rudder servo can be calculated from the step response of up to 100 seconds. Ans.: K = 20.5537, min{Jo} = 0.672

Figure 2.31. A simple yaw damper system (Prob. 2.8).

2.9. For the optimal problem described by

x1  x2 x2  u 5

,

x1 (0)  0

,

x2 (0)  1

min{ ( x12  5 x22  0.25u 2 )dt} 0

Controller Design Based on Parameter Optimization

83

i. Apply the method described in section 2.7 to show that it is not possible to reach the terminal states x1(5) = 1, and x2(5) = 1. ii. Modify the given optimal criterion to 5

min{[1000( x1 (t )  x1 (5)) 2  1000( x2 (t )  x2 (5)) 2 ] |t t   ( x12  5x22  0.25u 2 )dt} f

0

Show that the terminal states are reached within 10-3 accuracy. iii. Obtain the transfer function of the optimal controller for ii) above, and develop a Simulink model to simulate a closed-loop system. Ans.: J* = 4.6193. Figure 2.32 shows the optimal control and trajectories 2.10 A dynamic system is described by the nonlinear state equation [17]

x   0.2x  10 tanh( u), x(0)  5 Use the method of the vector control parameterization to determine u*(t), which minimizes the performance index tf

J (u )  10x 2 (t f )   (10x 2  u 2 )dt, t f  0.5 0

Find the corresponding controller transfer function and hence simulate the problem to determine the optimal trajectory.

Figure 2.32. Optimal control and trajectories (Prob. 2.9).

84

Optimal Control Engineering with MATLAB

2.11 Use genetic algorithm technique to solve the problem of Example 2.5. Show that the minimum value of ISE index is very near to that based on theoretical calculation, which is used in Example 2.5. Ans.: One possible transfer function is Gs ( s) 

( s  1.1742)( s  0.6117) ( s  3.5968)( s  1.512)

2.12 For the unstable plant described by the transfer function

Use genetic algorithm technique to determine the PID stabilizer, which minimizes the ITSE performance index. Hint: Use 1000 generations

LIST OF MATLAB PROGRAMS List 2.1 A Matlab program (Example 2.3) % The controller parameters are K and T clear all global T K t x v = [.1 1]; % initial estimated values for T and K [A, j] = fminsearch (@ two_par,v) optimset('TolX',1e-6); % the optimized values of K and T T = A(1); K=A(2); num = 500*[K*T K]; den = [conv([1 0], conv([1 16 100],[1 5]))] P = tf(num,den); G = feedback(P,1); % closed-loop with controller s = conv([1 0], conv([1 16 100],[1 5])); g = tf([500], [s]) gc = feedback(g,1); % closed-loop with negative unity feedback step (G, gc) function Jo = two_par(v) global T K t x T = v(1); K = v(2); B0 = 25e4*K^2; B1 = 25e4*K^2*T^2; b0 = 500*K; b1 = 500*K*T; a4 = 500*K; md = [500*K -180 1 0; 0 500*(1+K*T) -21 0; 0 -500*K 180 -1; 0 0 -500*(1+K*T) 21]; md0 = [500*(K*T+1) -180 1 0; 500*K 500*(1+K*T) -21 0; 0 -500*K 180 -1; 0 0 -500*(1+K*T) 21]; md1 = [500*K 500*(K*T+1) 1 0; 0 500*K -21 0; 0 0 180 -1; 0 0 -500*(1+K*T) 21]; d = det (md);d0 = det (md0); d1 = det (md1); Jo = (B0*d0+B1*d1-2*b1*b0*d)/(2*d*a4^2);

md = [500*K -180 1 0; 0 500*(1+K*T) -21 0; 0 -500*K 180 -1; 0 0 -500*(1+K*T) 21]; md0 = [500*(K*T+1) -180 1 0; 500*K Controller 500*(1+K*T) -21Based 0; Design on Parameter Optimization 0 -500*K 180 -1; 0 0 -500*(1+K*T) 21]; md1 = [500*K 500*(K*T+1) 1 0; 0 500*K -21 0; 0 0 180 -1; 0 0 -500*(1+K*T) 21]; d = det (md);d0 = det (md0); d1 = det (md1); Jo = (B0*d0+B1*d1-2*b1*b0*d)/(2*d*a4^2);

List 2.2 A Matlab program (Example 2.4) % Main program % The controller parameters are K and T % The index of performance is to minimize the ISE index for step input clear all global T K t x v = [ .1 1]; % initial estimated values for T and K [A,j] = fminsearch(@two_KT,v) optimset('TolX',1e-6); T = A(1); K = A(2); % the optimized values of K and T num = 500*[K*T K]; den = [conv ([1 0], conv ([1 16 100],[1 5]))] P = tf(num,den); G = feedback(P,1); % closed-loop with controller s = conv([1 0], conv([1 16 100],[1 5])); g = tf([500],[s]) gc = feedback(g,1); % closed-loop with negative unity feedback step(G, gc) function Jo=two_KT(v) global T K t x T = v(1); K = v(2); co = 500*(K*T+1); dHa = det([21 co 0 0; 1 180 500*K 0; 0 21 co 0; 0 1 180 500*K]); dHb = det([1 -81 11400 -250000; 1 180 500*K 0; 0 21 co 0; 0 1 180 500*K]); Jo = 0.5*dHb/dHa;

85

86

Optimal Control Engineering with MATLAB

List 2.3a A Matlab program (Example 2.5) % The controller parameters are z1, z2, p1, and p2 % The index of performance is to minimize ISE for unit-step clear all v = [0.5 0.5 2.0 0.5]; [P, j] = fminsearch(@Four_P,v); disp(' The Optimized Parameters') P disp(' The Minimal Index Value') j z1 = P(1); z2 = P(2); p1 = P(3); p2 = P(4); numc = conv([1 z1], [1 z2]); denc = conv([1 p1],[1 p2]); nums = [40]; dens = conv([1 0], conv([1 1],[1 4])); num = conv(numc, nums); den=conv(denc, dens); g = tf(num, den); gc = feedback(g,1); numr = [40 24 3.2]; denr = [1 9.02 24.18 16.48 0.32 0]; gr = tf(numr,denr); grc=feedback(gr,1); step(gc,'r',grc,'b'); system1 = allmargin(gc) system2 = allmargin(grc) function Jo = Four_P(v) z1 = v(1); z2 = v(2); p1 = v(3); p2= v(4); a1 = 40*z1*z2; a2=4*p1*p2+40*z1+40*z2; a3 = 4*p1+4*p2+5*p1*p2+40; a4 = 4+5*p1+5*p2+p1*p2; a5 = 5+p1+p2; a6 = 1;

Controller Design Based on Parameter Optimization

b1 = 4*p1*p2; b2 = 4*p1+4*p2+5*p1*p2; b3 = 4+5*p1+5*p2+p1*p2; b4 = 5+p1+p2; b5 = 1; A = [a1 a2 a3 a4 a5 a6]; B = [b1 b2 b3 b4 b5]; n = length(A) ; m = length(B); for i = 1: n-1 for j = 1: n-1 if ((2*j) = i Ha(i,j) = A(n-2*j+i); else Ha(i,j) = 0; end end end BR(1) = B(m)^2; BR(m) = (-1)^(m-1)*B(1)^2; for i = 2: m-1 s = 0; for j = 1: i-1 k1 = m+j-i+1; k2 = m-j-i+1; if k10 s = s + (-1)^j*(2*B(k1)*B(k2)); else end end BR(i) = (-1)^(i-1)*(B(m-i+1)^2+s); end Hb = [BR;Ha(2:n-1,1:n-1)]; dHa = det(Ha); dHb = det(Hb); Jo = 0.5*dHb/dHa;

87

88

Optimal Control Engineering with MATLAB

List 2.3b A Matlab program (Example 2.5) function Jo = Four_PC(v) z1 = v(1); z2 = v(2); p1 = v(3); p2 = v(4); a1 = 40*z1*z2; a2 = 4*p1*p2+40*z1+40*z2; a3 = 4*p1+4*p2+5*p1*p2+40; a4 = 4+5*p1+5*p2+p1*p2; a5 = 5+p1+p2; a6 = 1; b1 = 4*p1*p2; b2 = 4*p1+4*p2+5*p1*p2; b3 = 4+5*p1+5*p2+p1*p2; b4 = 5+p1+p2; b5 = 1; A = [a1 a2 a3 a4 a5 a6]; B = [b1 b2 b3 b4 b5]; n = length(A); m = length(B); num = [b5 b4 b3 b2 b1]; den =[ a6 a5 a4 a3 a2 a1]; E = tf(conv (num,[1 0]), den); tt = 0:.02:20; [e t] = step(E, tt); m = length(t); integd = (e.^2); h = t(2) - t(1); k = m-1; s = 0; s1 = 0; for i = 2:2:k s = s + integd(i); end k = k-1; for i = 3:2:k s1 = s1 + integd(i); end Jo = h*(integd(1)+integd(m)+4*s+2*s1)/3;

Controller Design Based on Parameter Optimization

List 2.4 A Matlab program (Example 2.8) clear all v = [10 10 10];% initial guess [P,j] = fminsearch(@PID_d,v); disp(' The Optimized Parameters') P disp(' The Minimal Index Value') j k1 = P(1); k2 = P(2); k3 = P(3); numc = [k1 k2 k3]; denc = [1 0]; nums = [1]; dens = [1 3.6 9]; numcs = conv(numc,nums); dencs = conv(denc,dens); g = tf(numcs,dencs); gc = feedback(g,1) % optimal control system Q = [5 10 15]; % three suboptimal control systems for i= 1:3 q = Q(1,i); numcr = [k1 k2 k3]/q; dencr = [1 0]; numcrs = conv(numcr,nums); dencrs = conv(dencr,dens); gr = tf(numcrs,dencrs); gcr = feedback(gr,1) tt = 0:.01:1; step(gc,gcr) hold on end function Jo = PID_d(v) k1 = v(1); k2 = v(2); k3 = v(3); a1 = k3; a2 = k2+9; a3 = k1+3.6; a4 = 1;

89

90

Optimal Control Engineering with MATLAB

b1 = 9; b2 = 3.6; b3 = 1; A = [a1 a2 a3 a4]; den = [a4 a3 a2 a1]; num = [b3 b2 b1 0]; E = tf(num,den); tt = 0:.01:5; [e t] = step(E,tt); h = t(2)-t(1); m = length(t); % Numerical computation of error derivative M = [-25 48 -36 16 -3; -3 -10 18 -6 1;1 -8 0 8 -1;-1 6 -18 10 3;3 -16 36 -48 25]; for i = 1:5:m-1 ec = [e(i);e(i+1);e(i+2);e(i+3);e(i+4)]; ed(i:i+4) = (M*ec)/(12*h); end ed(m) = (3*e(m-4) - 16*e(m-3) + 36*e(m-2) - 48*e(m-1) + 25*e(m))/(12*h); % Numerical computation of the index integd = t .*(abs(e)+abs(ed')); k = m-1; s = 0; s1 = 0; for i = 2:2:k s= s + integd(i); end k = k-1; for i = 3:2:k s1 = s1 + integd(i); end Jo = h*(integd(1)+integd(m)+4*s+2*s1)/3; % The stability test rce = roots(A); r=max(real(rce)); if r>=0|abs(Jo) 0, and k22 > 0, should be valid. The algebraic system is given by

k12  4  0 2

k 2 k12  0 k11  2k12  k12 k 22  0 2 2k12  4k 22  k 22 0

2k 2  k 2 k 22  k1  0 k 22  0 The solution is

k  k1  k 2  0 k122  4 k 222  4k 22  2k12  0

or,

k 22  [4  16  8k12 ] / 2  2( 2  1) k11  4 2

286

Optimal Control Engineering with MATLAB Hence, the optimal control and V* are given as

u * (t )  2 x1 (t )  2( 2  1) x2 (t )

V * ( x(t ))  4 2 x12 (t )  4 x1 (t ) x2 (t )  2( 2  1) x22 (t )  0 The above solution is adopted because it satisfies the positiveness condition, i.e.

4 2 2  det(   )  0.68629  0  2 2( 2  1) It can be easily shown that at the steady-state the closed-loop optimal control system is stable with a natural frequency of a value, wn  2 and a damping coefficient of a value,   1 . Furthermore, with initial values x1(0) = 1, and x2(0) = 0, the optimal trajectories x1*(t) and x2*(t) can be obtained as

x1 * (t )  e 

2t



2te 

x2 * (t )  2te 

2t

2t

In fact, this steady-state solution is valid for t ≥ 3.5 as shown in the Figure 6.12, where it shows the result of numerical backward integration of the differential system with step of integration equal to h = 0.0001(see the MATLAB program given in List 6.3). Figure 6.12a shows very near values of parameters k‟s to those obtained by the steady-state solution. Furthermore, these results indicate the correctness of the used numerical calculation, which can be used for fixed time problems. If, for example, tf = 5 units of time, then the steady-state solution is accepted, and the optimal control and trajectories are as shown in Figure 6.12b. The optimum value of the optimal criterion is approximately 5.6. For the sake of completeness let us assume that the given problem is of finite time one, say for tf = 0.5 (it should be noted that this is not a case of a free terminal time). In this case, the previous numerical calculation can be carried out to get the results shown in Figure 6.13. A time-varying optimal control structure can be now attained by fitting the numerical data to time functions. In our example, the following functions can be established as polynomials of third-order.

k11 (t )  0.2918t 3  0.324t 2  1.2193t  1.9852 k12 (t )  0.9067t 3  0.5288t 2  1.2193t  0.3639 k 22 (t )  0.1929t 3  0.583t 2  0.4091t  0.0836

Dynamic Programming

287

Figure 6.12. a The k's parameters of the parametric expansion (Example 6.8).

Figure 6.12. b Optimal control and trajectories (Example 6.8).

The values of k2, which affects the optimal control, are very small and one can but k2 = 0 as in the steady-state solution. In this case, the optimum value of the optimal criterion is approximately 2. If the time-invariant optimal controller, which is obtained with the steadystate case, is used for tf = 0.5 then the optimum value is 3.1, and hence less minimum than with time-varying optimal controller. Therefore, the time-invariant controller can be accepted

288

Optimal Control Engineering with MATLAB

as an optimal one when final time is less than a certain value; practically, it may be accepted as a suboptimal solution. As a final note, the reader can use the program in List 6.4 (and the Simulink model shown of the time-varying controller in Figure 6.14) to check the answers, we got in this example.

Figure 6.13. The k's parameters for the interval [0, 0.5] (Example 6.8).

Figure 6.14. Simulink model of a time-varying controller (Example 6.8).

Dynamic Programming

289

6.6. CROSS-PRODUCT TERMS Based on the benefits we got from previous discussions and examples, we would like to introduce again the linear-quadratic control problem in the most general form or what is called "the cross-product terms problem"; see example 3.5. In many practical cases of nonlinear and/or non-quadratic problems, the linearization and quadraticization of the optimal criterion about the optimal trajectory yield a linear time-varying system with a quadratic criterion containing such cross-product terms [16,17]. Let us consider a linear time-varying MIMO system and a quadratic optimal criterion with finite terminal time described as

x (t )  A(t ) x(t )  B (t ) u ,

x  n ,

u  m

V  min{x T (t f ) M (t f ) x(t f )

(6.36)

u

tf

  [ x T Q x  u T R u  2u T N x  2Tx  2 Su] dt} t0

where Q is positive semi-definite and R is symmetric positive definite. If it is assumed that the pair (A, B) is stabilizable, then the parametric expansion is guessed as a possible form of the solution for the optimal criterion. We put it in the form

V ( x(t ))  K 0 (t )  2 xT K1 (t )  xT K 2 (t ) x

(6.37)

where K2 is asserted to be symmetric and the following positive semi-definite constraint holding

Q  NR 1 N T  0

(6.38)

Following the usual derivation of solving the HJB equation with some vector calculus manipulations, we obtain the optimal control vector as

u * (t )   R 1[( N  B T K 2 ) x  ( S T  B T K1 )]

(6.39)

where the matrices K1 and K2 are obtained from solving the following matrix first-order differential equations.

K1(t )  ( N  BT K 2 )T R 1 ( S  BT K1 )  T T  AT K1 K 2 (t )  ( N  BT K 2 )T R 1 ( N  BT K 2 )  Q  K 2 A  AT K 2 with terminal conditions K1(tf) = 0, and K2(tf) = M(tf).

(6.40) (6.41)

290

Optimal Control Engineering with MATLAB

In applications, it is found that most cases include only cross-product terms between the state vector x and the input vector u, i.e. only N is a nonzero matrix while matrices T and S do not exist (in this case K0 and K1 are of no meaning to assume in the V function of Equation (6.37)). Moreover, if the terminal time is infinite and the steady-state solution is required, the optimal control will be given by

u * (t )   R 1 ( N  BT K 2 ) x (t )

(6.42)

( N  B T K 2 ) T R 1 ( N  B T K 2 )  Q  K 2 A  AT K 2  0

(6.43)

The solution obtained from Equations (6.42) and (6.43) is exactly what the MATLAB function lqr does; an error message is also reported whenever the necessary condition and the constraint are not satisfied. In fact, it represents a special solution of the HJB equation.

6.7. DIFFERENTIAL DYNAMIC PROGRAMMING DDP In previous sections, the application of dynamic programming considers both discrete and continuous systems of relatively law-order. Although, a global optimal solution of several problems has been reached, the challenges of the curse of dimensionality with discrete-time systems and the unavailability of a general methodology for solving the HJB equation with continuous systems preclude the use of the conventional dynamic programming for highorders and/or nonlinear systems. D. Jacobson and D. Mayne [18] were the pioneers of one of the very interesting local methods that solve a big class of optimal problems, the differential dynamic programming DDP. The DDP is an iterative procedure based on local expansion to a second-order of the optimal cost function. The method is local (whenever the convexity conditions are satisfied), in the sense that it maintains a representation of a single trajectory and improves it locally. However, the improvement itself is based on dynamic programming – within a "tube" around the current trajectory. In fact, an ideal blend of the advantages of local and global methods is provided by DDP. The method can be applied to discrete-time or continuous-time optimal control problems. In the last decade, DDP was applied to the hard problem of motor control and robot locomotion, with remarkable success [19, 20]. In what follows, we will consider first a discrete formulation of DDP. Let us consider the discrete unconstrained free end-point problem defined as

x(k  1)  f ( x(k ), u(k ), k )

(6.44)

N 1

min{V ( x(0), u )   ( x( N ))   L( x(k ), u (k ), k )} k 0

where x(0) is given, and the control time is discretized to N steps.

(6.45)

Dynamic Programming

291

~(k ) approximates the optimum control and ~ Let a nominal control u x (k ) be the ~(k ) in the corresponding approximation of the optimum trajectory obtained by substituting u dynamic Equation (6.44). To improve the approximation, the DDP method suggests a second~(k ) and ~ order expansions of L and f around the nominal estimates u x (k ) , (for simple writing we drop the argument k from the partial derivatives f x , f xx , Lx , Lxx that are calculated at

u~(k ) and ~ x (k )) , thus L( ~ x (k )   x(k ), u~(k )   u (k ))  L( ~ x (k ), u~(k ))   x T (k ) Lx   u T (k ) Lu  1  x   Lxx 2  u   Lux T



Lxu   x   O3 Luu   u 

1  x(k )  f xx  x(k  1)   x (k ) f x   u (k ) f u   2  u (k )  f ux T

T

T

(6.46)

f xu   x(k )  O3    f uu   u (k ) (6.47)

where  x(k )  x(k )  ~ x (k ) , and the cubic terms O 3 ( x,  u ) will be neglected. Similarly, a second-order expansion of the unknown cost function V ( x(k )) about the nominal trajectory can be carried out

1 V ( x(k ))  V ( ~ x (k )   x(k ))  V ( ~ x (k ))   x T (k )Vv   x T (k )Vxx  x(k ) (6.48) 2 where the partial derivatives Vx , Vxx are calculated on  ~ x (k ) . Let us define the variation of the cost function by

 V (~ x (k ), u~(k ))  V (~ x (k ), u~)  V (~ x (k )) , Then we can rewrite Equation (6.48) as

1 V ( x(k ))  V ( ~ x (k ), u~)   V ( ~ x (k ), u~(k ))  Vx x(k )  Vxx  x(k ) 2

(6.49)

The application of the dynamic programming relation (Equation (6.1)) for the cost function L, and the trajectory variation x given in Equations (6.46) and (6.47), results in

Optimal Control Engineering with MATLAB

292

1 V (~ x (k ), u~)  V ( ~ x (k ), u~(k ))  V x x(k )  V xx  x(k )  min{L( ~ x (k ), u~ (k )) u ( k ) 2  Lx  1  x(k )  Lxx  L   x(k )  2  u (k )  L    ux  u T

Lxu   x(k )  V (~ x (k  1), u~)   V ( ~ x (k  1), u~ )    Luu   u (k )

1  V x (k  1) x(k  1)   x T (k  1)V xx (k  1) x(k  1)} 2 (6.50) To minimize equation 6.50 based on the principle of dynamic programming; the definition of Hamiltonian function H ( x(k ), u(k ),V (k  1) could be used to simplify Equation (6.50). In addition, on the assumption that all partial derivatives are calculated on ~ x (k ), u~(k ) , the DDP algorithm for improving the nominal control toward the optimal one could be stated as

 u(k )  C 1 (k )[H u (k )  B(k ) x(k )] where

(6.51)

~ x (k ) is given by Equation (6.47), and C (k )  H uu (k )  f x (k )Vxx (k  1) f x (k ) T

B (k )  H ux (k )  f u (k )Vxx (k  1) f x (k ) T

Vxx (k )  A(k )  B T (k )C 1 (k ) B(k ) A(k )  H xx (k )  f x (k )Vxx (k  1) f x (k ) T

(6.52a)

Vx (k )  H x (k )  B T (k )C 1 (k ) H u (k ) H x (k )  [ Lx (k )  Vx (k  1) f x (k )]T T

T

H u (k )  [ Lu (k )  Vx (k  1) f u (k )]T T

T

The functions of k are conveniently calculated by stepping backwards from the boundary values

Vx ( N )   x ( ~ x ( N )) V ( N )   (~ x ( N )) xx

(6.52b)

xx

To permit the possibility of restricting the magnitude of the iterative corrections of the control, Equation (6.50) has to be slightly changed to

 u(k )  C 1 (k )[ H u (k )  B(k ) x(k )], 0    1

(6.53)

Dynamic Programming

293

Jacobson [21] showed that the evaluation of the partial derivatives could be carried out ~(k )  x * (k ) , where the control variation x * (k ) is along a trajectory computed by u given by a modified version of Equation (6.53) as

 u * (k )  C 1 (k )[ H u (k )], 0    1

(6.54)

Note that the computation with Equation (6.54) is simpler than with (6.53). However, during the iterative procedure if, it happens that ε = 1 results in a higher cost value than in the previous iteration, then the bulk of computation per step has to be repeated with smaller value of ε. The convexity requirement is satisfied as indicated in Equations (6.52), if both Huu(k) and Vxx(k+1) are positive semi-definite and one of them at least is positive definite. For discretetime system, the DDP iterative computation sequences can be summarized as follows: 1. Assume a nominal approximated solution ũ(k), k = 0,1…N -1; usually a (suitable) constant value for the whole control interval. Also set ε ≤ 1. 2. Compute the nominal trajectory (k) by means of the given dynamic (Equation (6.44)) and the nominal cost V( , ũ) from Equation (6.53) (or 6.54). 3. Computing V x and V xx by stepping backward Equations (6.51) starting with terminal

4. 5. 6. 7.

conditions given in Equations (6.52); checking whether the convexity conditions have been satisfied otherwise the nominal control has to be changed. Compute δu*(k) by Equation (6.54) and apply the correction to the nominal control, i.e., u(k) = ũ(k) + δu*(k) as the next approximation in step number 2. Compute the new trajectory and the new cost. If the old cost is greater than the new cost then set ε = ε /2 and return to 2 above but with the latest u(k), else Check the stopping criterion or end the prescribed number of iterations; the stopping criterion could be simply a prescribed accuracy of the minimization process.

Example 6.9. Consider the first-order nonlinear discrete-time system

x(k  1)  (1  0.02x(k ))x(k )  0.01u(k ) It is required to obtain the optimal control sequence, which minimizes the cost function

0.5 99 2 V  0.2 x 2 (100)  ( )  x (k )  u 2 (k ) 101 k 0 The partial derivatives (as functions of k) of the dynamic system and of the Hamiltonian function with respect to the state and control are

294

Optimal Control Engineering with MATLAB

f x (k )  1  0.04 x(k ) f u (k )  0.01 0.5 x( k )  (1  0.04 x(k ))Vx (k  1) 101 0.5 H x x (k )  2  0.04Vx (k  1) 101 0.5 H u (k )  2 u (k )  0.01Vx (k  1) 101 0.5 H uu (k )  2 101 H x (k )  2

The backward calculation starts with

V x (100)  0.4 x(100)

V xx (100)  0.4

~(k )  0.1, k  0,1,2,,99 . Following the The nominal control is selected to be u given DDP algorithm, a MATLAB program is written as in List 6.5. When the value of ε is set initially to zero, after 10 iterations the discrete optimal control and trajectory, and the cost minimization process are as depicted in Figures 6.15 and 6.16.

Figure 6.15. Optimal control and trajectory (Example 6.9).

The initial and final optimal control values are respectively -0.229 and -0.102 and the state reaches 0.253 at the terminal time. The cost function decreases sharply, and it terminates to 0.17044 with an accuracy of 10-8 .Note also that the value of ε is not reduced to smaller values; the case may be different with other parameters.

Dynamic Programming

295

Figure 6.16. Cost function versus iterations (Example 6.9).

The reader could solve with grater terminal k say 149 (in the program N = 151) and it will be found that ε has to be reduced to 0.125 and that the state will reach the value of 0.158. However, as the terminal k increases the design parameters (number of iterations and the initial value of ε) have to be changed. For example, for N = 201, an accepted optimal response will be obtained with 40 iterations and initial ε equal to 0.25. The trajectory reaches 0.0847 and the minimum cost becomes equal to 0.09. Finally, the required number of iterations depends on the nominal value of control; as the nominal value is far from the suitable value, more iterations will be required. Incorrect nominal values of control lead to the breaking of the convexity conditions; for example, the ~(k )  0.4 . convexity conditions are not satisfied for N = 201, and u For the continuous optimal control problem

x (t )  f ( x, u , t ), V  min{ ( x (t f )) 

x ( 0)  x 0 tf

t L( x, u, t ) dt}

(6.55)

0

Using the limiting process (the step of discretization h tends to zero) a corresponding set of the DDP equations were derived by Mayne [22], as

296

Optimal Control Engineering with MATLAB

H ( x, u , V x , t )  L ( x, u , t )  V x f ( x, u , t ) T

 u (t )   H uu 1 [ H u  ( H ux  f u T V xx ) x(t )] dV 1 T  x  H x  H uu [ H ux  f u V xx ]T H u dt dV T  xx  H xx  f x V xx  V xx f x  dt T 1 T  [ H ux  f u V xx ]T H uu [ H ux  f u V xx ] V x (t f )   x ( x(t f ))

(6.56)

V xx (t f )   xx ( x(t f )) u  u~   u Although, the last Equation of (6.56) computes the control by successive approximation starting with nominal control, the DDP algorithm is originally derived such as to compute the control using

u(t )  u~(t )   u * (t )   (t ) x(t )  u * (t )   (t )  x(t )

(6.57)

where

   H uu 1 ( H ux  f u TVxx )

(6.58)

Clearly, the two coupled ordinary differential equations should be integrated in backward procedure starting with the given boundary conditions and suitable nominal control (usually constant value over the whole control interval). Obviously, most nonlinear problems require numerical integration. This could be carried out by defining a continuous control via applying continuous-time integration such as RungaKutta method or Euler method with sufficient small integration step. However, it is necessary that the solution of Vxx Ricatti equation should be bounded to obtain bounded Vx solution and in turn a solution convergence. The sufficient conditions that secure boundedness of Ricatti equations are that 1

H xx  H ux H uu H ux  0 T

1

H uu  0

 xx  0

(6.59)

Dynamic Programming

297

Some of the DDP characteristics are [18] 1. It exhibits one-step convergence on linear-quadratic cost problems; it very important results for linearize systems. 2. In a neighborhood of the optimum, convergence is rapid for non-linear-quadratic cost problems. 1

x , u,Vx ; t ) at the minimizing u*, and 3. It requires only a positive definite H uu ( ~ therefore, the algorithm can handle a larger class of non-linear-quadratic cost problems. 4. It is capable of computing the optimal control even with problems that have an unbounded solution of the Ricatti equations along some nominal trajectories. The next example considers a continuously stirred tank reactor CSTR [23], described by second-order nonlinear dynamic.

Example 6.10. For the nonlinear CSTR

 x    0.01x1  0.338 0.02   x1  0.02  1       u,  x2   0.05x1  0.159  0.03  x2   0 

 1.5 x(0)     3 

It is required to find the optimal control, which regulates (to bring the trajectories to zero with zero final optimal control) the reactor while minimizing the optimal criterion T

min{ J   ( x T Qx  Ru 2 ) dt}, 0

10 0 Q  , R  1, t  [0,300]  0 1 In such simple nonlinear-quadratic problem, the long control time replaces the infinite horizon time (the steady-state case). The Hamiltonian function and its partial derivatives are

298

Optimal Control Engineering with MATLAB

H  10x1  x2  u 2  Vx1 (0.01x1  0.338x1  0.02x2  0.02u ) 2

2

2

 Vx2 (0.05x1  0.159x1  0.03x2 ) 2

20x1  Vx1 (0.02x1  0.338)  Vx2 (0.1x1  0.159) Hx    2 x2  0.02Vx1  0.03Vx2   20  0.02Vx1  0.1Vx2 0 H xx   0 2  H u  2u  0.02Vx1

 u*  0.01Vx1

H uu  2  0 H ux  0 The dynamic partial derivatives are

 0.02x1  0.338 0.02  fx   ,  0.1x1  0.159  0.03

0.02 fu     0 

The backward calculation starts with

0  Vx (t f )    0  0 0  Vxx (t f )    0 0  The sufficient conditions for a bounded solution have to be checked first. The second and third conditions are always satisfied and since Hux is zero the first condition will be satisfied if and only if Hxx is positive semi-definite, or

20  0.02Vx1  0.1Vx2  0 It turns out that unbounded solution may take place. The nominal control is chosen to be u~(t )  0, t [0,300] . The two differential equations are integrated by Euler method with a step of integration h = 0.1. A MATLAB program (shown in List 6.6) is written to execute the continuous DDP algorithm defined by Equations (6.56). The solution converges very rapidly. Figure 6.17 shows "the optimal control tube" and Figure 6.18 shows the final results (after 20 iterations) of the optimal control and trajectories. Actually, with zero nominal control, the final results can be obtained in a less number of iterations; with different nominal values say u~(t )  1, the final results would not be different but the tube has a different picture as depicted in Figure 6.19. Furthermore, as shown, the CSTR is regulated with a minimal cost value of 169.541 within 10-4 accuracy.

Dynamic Programming

299

Figure 6.17. Optimal control tube (Example 6.10).

Figure 6.18. Optimal control and trajectories (Example 6.10).

The reader may run the program while computing the improvements in control according to Equations (6.57-58) to get the same optimal control and trajectories and the same minimal cost value. Moreover, he can try obtaining the same results for different nominal values say in the range -3 ≤ u0 ≤ 3 and noting the increment of the required number of iterations (for u0 =3, 30 iterations are required). Finally, he may wish to check for an unbounded solution by

300

Optimal Control Engineering with MATLAB

setting the nominal control u0 equal to say 4 (or by decreasing the control weight value R say to 0.5). Finally, a suboptimal state feedback controller can be found (see Example 4.12) as

u * (t )  0.5123x1 (t )  0.2201x2 (t ) Furthermore, it is of interest to implement the suboptimal controller via the setting of a Simulink model (see Figure 4.20). Obviously, this implementation will result slightly higher cost value than the optimal value.

For problems with terminal constraints, bang-bang control problems and other related issues including comparison with similar but different approaches, the reader could consult the book of Jacobson and Mayne [18]. A smart technique (named step size adjustment method) for solving the unboundedness problem is also suggested in the reference.

Figure 6.19. The optimal control tube with

u~(t )  1 (Example 6.10).

PROBLEMS 6.1. Use the dynamic programming approach to solve the following constrained optimal control problem

Dynamic Programming

301

x   u (t ), x(0)  x0 , x(5)  2, | u (t ) | 1 5

J   (2 x 2  u 2 ) dt 0

Hint: solve as it is illustrated in Example 6.1 6.2. Consider the first-order discrete-time system

x(k  1)  ax(k )  bu(k ) It is required to minimize the optimal criterion

V 

k 2

x

2

(k )  ru 2 (k ) , r > 0

k 0

Apply the dynamic programming approach to show that

ra 2 ) r  b2 u * (0)   x ( 0) ra 2 2 r  b (1  ) r  b2 ab(1 

What is the expression of u*(1)? 6.3 Consider a scalar discrete-time system defined by

x(k  1)  2x(k )  3u(k ),

x(0)  4

Apply the principle of optimality to determine the control and state sequences, which minimize the optimal criterion

J  ( x(2)  10) 2 

1 1  [ x 2 (k )  u 2 (k )] 2 k 0

Answer: u*(0) = 2.0153, u*(1) = -1.9237 6.4 A second-order discrete-time system is described by

0.0785  x1 (k )  0.00308  x1 (k  1)   0.997  x (k  1)   0.0785 0.997   x (k )   0.0785  u (k )  2   2    

302

Optimal Control Engineering with MATLAB

Determine the discrete optimal feedback control u*(k), which minimizes the optimal criterion 

V   [ x12 (k )  x22 (k )  u 2 (k )] k 0

Answer: u * (k )  0.34x1 (k )  1.29x2 (k ) 6.5 For the voltage regulator problem [16], the following system matrices and optimal criterion are found

 0.2  0  A 0   0  0 

C  1 0

0.5

0

0

 0.5 0

1.6 100  7 0

0 600 7  25

0

0

0

0

0

0

0,

0  0  0 ,  75   10

0 0   B   0 ,   0 30



min{J   [ x12 (t )  u 2 (t )] dt} 0

In the reference, the solution of Ricatti equation gives the following optimal gain vector

k  0.9243 0.1711 0.0161 0.0492 0.2644 Modify the MATLAB program given in List 6.2 to obtain the optimal gain vector D (see Equation (6.14)). Find the suitable simulation time (the N parameter in the program) to reach the result of k given above. 6.6 For the first-order system defined by

x(t )  3x(t )  u(t ),

x(0)  1

It is required to transfer the system to x(T) = 0.5 while minimizing the functional T

V 

1 [8 x 2 (t )  u 2 (t )] dt 2 0

a) Show that the solution of the corresponding HJB equation yields to the following Ricatti equation

Dynamic Programming

P   6 P  2 P 2  4  0,

303

P (T )  0

b) Determine u*(x, t). c) Simulate the optimal solution to calculate T which satisfies the terminal condition. 6.7. Consider the following numerical optimal control problem; a first-order nonlinear system with non-quadratic optimal criterion

x(t )  8 x(t )  2u (t ),

x(0)  3



V  min{ [3x 3 (t )  0.1u 2 (t )] dt} 0

Determine the optimal control u*(t) and simulate to determine V*. Hint: Use Equation (6.31) 6.8. Consider the scalar optimal control problem

x1 (t )   x(t )  u (t )  1,

x(0)  1

1

V (u )  x 2 (1)   u 2 (t ) dt 0

Show that the optimal feedback control is given by

u * (t )  

4 (e 4  e 2(t 1) ) 4e 4  x(t ) 7e 4t  3e 4 7e 4t  3e 4

Use the parametric expansion

V * (t , x)  k (t )  k1 (t ) x(t ) 

1 k 2 (t ) x 2 (t ) . 2

6.9. Consider the scalar nonlinear optimal control problem

x1 (t )  x(t )u (t ),

x(0)  1

1

V (u )  x 2 (1)   [ x(t )u (t ) ]2 dt 0

Use V (t, x) = p (t) x2(t) to determine the optimal control laws u*(t) and u*(x). Simulate to verify your solution. Answer: x*(1) = 0.5, V*=0.5.

Optimal Control Engineering with MATLAB

304

6.10. For double integral system defined by the state equation

x1 (t )  x2 (t ) x2 (t )  u (t ) Solve the HJB equation such as to determine the optimal state feedback control, which minimizes the cost function 

1 V   [ x12 (t )  2bx1 (t ) x2 (t )  ax22 (t )  u 2 (t )] dt, 20

(a  b 2 )  0

and makes the closed-loop system stable. Hint: Use V  Ax1 (t )  Bx1 (t ) x2 (t )  Cx 2 (t ) 2

2

6.11. For the optimal control problem

x1 (t )  x 2 (t ) x 2 (t )  ax2 (t )  bu(t ) 

V  min{ [ 1 x12 (t )   2 x 22 (t )   3 u 2 (t )] dt } 0

a) To determine the optimal control u*, show that it is always sufficient to admit

V *  k11 x12 (t )  2k12 x1 (t ) x2 (t )  k 22 x22 (t ) b) Determine the optimal control. 6.12. Consider the following optimal control problem

x1 (t )  x 2 (t ) x 2 (t )  x1 (t )  x 2 (t )  u (t ) tf

1 V  min{  [ x12 (t )  u 2 (t )] dt} 2 0 a) For infinite horizon case, determine the constant feedback control law, which makes the closed-loop system stable. b) If tf is equal to 1, determine the time-varying optimal control.

Dynamic Programming

305

Hint: solve as it is illustrated in Example 6.8. 6.13. Consider the following optimal control problem

x1 (t )  x 2 (t ) x 2 (t )   x 2 (t )  u (t ) V  min{



1 [ x12 (t )  x 22 (t )  u 2 (t )] dt}  20

If we assume that V *  k11 x1 (t )  2k12 x1 (t ) x2 (t )  k 22 x2 (t ) , show that the 2

2

solution of the HJB equations leads to three possible candidate solutions of the optimal control law. Determine the stabilizer one, and check whether it satisfies the boundary condition J*(x (∞)) = 0. 6.14. For the first-order linear system

x(t )  0.5x(t )  u(t ),

x(0)  1

Find the optimal control u* over t in the interval [0, 1] which minimize the performance index (of cross-product term) 1

V   (0.625x 2  0.5 xu  0.5u 2 ) dt 0

6.15. Apply DDP algorithm to solve the discrete optimal control problem [5]

x( k  1)  0.99 x( k )  0.5 tanh( u ( k )), V ( x(0), u ) 

x ( 0)  5

9

115 2 1 x (10)  [10 x 2 ( k )  u 2 ( k )]  11 22 k 0

6.16. Consider the nonlinear optimal control problem

x '   x 3  u ; x(0)  1 V  min[

1

1 ( x 2  u 2 )dt ]  20

Apply the discrete DDP to determine the optimal control and trajectory. 6.17. Apply continuous DDP version to determine the optimal control and trajectories.

306

Optimal Control Engineering with MATLAB

 x1  2 x1  x 22 u, x1 (0)  1  x 2  2 x13  x 2  2u, x 2 (0)  0 min{J 

0.5

0 ( x

4 1

 2 x 22  u 2 ) dt}

Hint: use zero nominal control and 0.0005 step of integration.

LIST OF MATLAB PROGRAMS List 6.1 Matlab program (Example 6.2) clear all a=1;b=3;c=0.25;x0=2;N=10;h=0.1; N=N+1; f=exp(-a*h);e=(-b/a)*(exp(-a*h)-1); dN=0;g(N)=0; % Starting the reverse iterative procedure of the dynamic programming for k=N-1:-1:1 m(k)=k; d(k)=e*f*g(k+1)/(g(k+1)*e^2+c*h); g(k)=(1+c*d(k)^2)*h+g(k+1)*(f-e*d(k))^2; end dk=flipud([d dN]'); gk=flipud(g'); mk=flipud([0 m]'); [mk dk gk] cp0=1; % Starting forward calculation of the optimal trajectory for i=1:N-2 cp(i)=cp0*(f-e*d(i)); cp0=cp(i); x(i)=cp(i)*x0; end xv=[x0 x];u=[-d.*xv]' [xv' u] T=0:1:N-1; tt=0:0.01:N-1; u=[u; 0]; y=0; % to plot the x-axis plot(tt,y,'k') hold on stairs(T,u,'k') hold on stairs(T,[xv 0])

Dynamic Programming

List 6.2 A Matlab program (Example 6.3) clear all a=[-0.37 .0123 .00055 -1;0 0 1 0;-6.37 0 -.23 .0618;1.25 0 .016 -.0457]; b=[.00084 .000236;0 0;.08 .804;-.0862 -.0665]; c=[0 1 0 0;0 0 1 0;0 0 0 1]; d=0; x0=[1;1;-1;-1];x00=x0; A=eye(4);H=eye(2); h=0.001; sys=ss(a,b,c,d); sysd=c2d(sys,h); [f,e]=ssdata(sysd) N=6001; G=zeros(4,4); for k=N-1:-1:1 m(k)=k; K=G*e*inv(H+e'*G*e); G=A+f'*(G-K*e'*G)*f; D=K'*f; x1=(f-e*D)*x0;% the optimal closed system x1v(k)=x1(1,:);x2v(k)=x1(2,:);x3v(k)=x1(3,:);x4v(k)=x1(4,:); x0=x1; end D % the initial gain matrix X=[x00';flipud([x1v' x2v' x3v' x4v'])]; t=h*[0 m];

plot(t, X(1:end,1).‟k‟, t, X(1:end,2),‟k‟, t, X(1:end,3),‟k‟,t,X(1:end,4),‟k‟) List 6.3 Matlab function (Example 6.5) function pt=dynamic(V) t=V(1);a=V(2);b=V(3);c1=V(4);c2=V(5);M=V(6);T=V(7); B=2*a*c2/(b^2);C=-c1*c2/(b^2); m1=0.5*(-B-sqrt(B^2-4*C));m2=0.5*(-B+sqrt(B^2-4*C)); k=(m1-M)/(m2*exp((m1-m2)*T*b^2/c2)); er=(m1-m2)*b^2/c2;

pt=-(m1-m2*k*exp(er*t)./(1-k*exp(er*t));

307

308

Optimal Control Engineering with MATLAB

List 6.4 A Matlab program (Example 6.8) clear all h=.0001;% step of integration T=3.5;% total time of integration M=T/h; % boundary conditions of the k-parameters k11(M)=0;k1(M)=0;k12(M)=0;k22(M)=0;k2(M)=0;k(M)=0; % Back integration of equations using Euler numerical method % for each parameter of the parametric expansion for i=M:-1:2 a=-.5*h+k12(i)^2-4; k11(i-1)=k11(i)-a*h; b=-0.5*h+k2(i)*k12(i); k1(i-1)=k1(i)-b*h; c1=k12(i)*k22(i)+2*k12(i)-k11(i); c2=-.5*h+(k12(i)-.5*c1*h)*k22(i)+2*(k12(i)-.5*c1*h)-k11(i); k12(i-1)=k12(i)-c2*h; d1=k22(i)^2+4*k22(i)-2*k12(i); d2=-.5*h+(k22(i)-.5*d1*h)^2+4*(k22(i)-.5*d1*h)-2*k12(i); k22(i-1)=k22(i)-d2*h; e1=k1(i)-(2+k22(i))*k2(i); e2=-.5*h+k1(i)-(2-k22(i))*(k2(i)-.5*e1*h); k2(i-1)=k2(i)-e2*h; f=-.5*h+k2(i)^2; k(i-1)=k(i)-f*h; end t=0:1:(M-1); tt=t*h; % plots of the k's parameters of the parametric expansion figure(1) plot(tt,k11,'k',tt,k1,'k',tt,k12,'k',tt,k22,'k',tt,k2,'k',tt,k,'k') xlabel('time') ylabel( 'k-parametric expansion') text(2.5,4,'k11'),text(1.7,1.9,'k12'),text(1.35,1,'k22'),text(0.6,.2,'k,k1,k2') %the optimal trajectories plots for steady-state case st=0:.01:5; x1=exp(-sqrt(2).*st)+sqrt(2).*st.*exp(-sqrt(2).*st); x2=-2*st.*exp(-sqrt(2).*st); uopt=-2*x1-2*(sqrt(2)-1)*x2;

Dynamic Programming

List 6.4 A Matlab program cont. figure (2) plot(st,uopt,'k',st,x1,'k',st,x2,'k') xlabel('time') ylabel( 'u*(t), x*1(t), x*2(t)') text(1,0.7,'x*1(t)'),text(1,-0.3,'x2*(t)'),text(1,-0.85,'u*2(t)') %fitting for k's parameters as a function of time p_k11=polyfit(tt,k11,3) p_k12=polyfit(tt,k12,3) p_k22=polyfit(tt,k22,3) p_k2=polyfit(tt,k2,3) List 6.5 A Matlab program (Example 6.9) % DDP for discrete nonlinear system clear all x0=1;S=0.2;qr=0.5; % problem parameters N=101; u0=-.1; % nominal control epslon=.5; u=u0*ones(N-1,1); NI=12; % number of iterations for m=1:NI % calculation of nominal trajectory for i=1:N-1 x(1)=x0; x(i+1,:)=(1-0.02*x(i,:))*x(i,:)+0.01*u(i,:); xn(i)=x(i); end Vx(N)=2*S*x(end); Vxx(N)=2*S; dx(N-1)=0; % cost function calculation pv=0.2*x(N)^2+(qr/N)*sum(x(1:N-1).^2+u(1:N-1).^2); if m==1 iv=pv; else,end % Backward Calculation for k=N-1:-1:1 fx(k)=1-0.04*x(k); fu(k)=0.01; hx(k)=2*(qr/N)*x(k)+(1-0.04*x(k))*Vx(k+1); hxx(k)=2*(qr/N)-0.04*Vx(k+1); hu(k)=2*(qr/N)*u(k)+0.01*Vx(k+1); huu(k)=2*qr/N; hux(k)=0; A(k)=hxx(k)+fx(k)'*Vxx(k+1)*fx(k); B(k)=hux(k)+fu(k)'*Vxx(k+1)*fx(k); C(k)=huu(k)+fu(k)'*Vxx(k+1)*fu(k); vx(k)=hx(k)+B(k)'*inv(C(k))*hu(k); Vxx(k)=A(k)+B(k)'*inv(C(k))*B(k);

309

310

Optimal Control Engineering with MATLAB

List 6.5 A Matlab program cont. % The convexity conditions if huu(k)