# Patent

## Title
Control system design method

## Abstract

A control system design method and concomitant control system comprising representing a physical apparatus to be controlled as a Hamiltonian system, determining elements of the Hamiltonian system representation which are power generators, power dissipators, and power storage devices, analyzing stability and performance of the Hamiltonian system based on the results of the determining step and determining necessary and sufficient conditions for stability of the Hamiltonian system, creating a stable control system based on the results of the analyzing step, and employing the resulting control system to control the physical apparatus.

## Background
CROSS-REFERENCE TO RELATED APPLICATIONS
This application claims priority to and the benefit of the filing of U.S. Provisional Patent Application Ser. No. 60/896,043, entitled “Exergy and Irreversible Entropy Production Thermodynamic Concepts for Control System Design”, filed on Mar. 21, 2007, and the specification and claims thereof are incorporated herein by reference.


STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT
The Government has rights to this invention pursuant to Contract No. DE-AC04-94AL85000 awarded by the U.S. Department of Energy.


INCORPORATION BY REFERENCE OF MATERIAL SUBMITTED ON A COMPACT DISC
Not Applicable.
COPYRIGHTED MATERIAL
Not Applicable.
BACKGROUND OF THE INVENTION
1. Field of the Invention (Technical Field)
The present invention relates to control system design methods.
2. Description of Related Art
Note that the following discussion refers to a number of publications by author(s) and year of publication, and that due to recent publication dates certain publications are not to be considered as prior art vis-a-vis the present invention. Discussion of such publications herein is given for more complete background and is not to be construed as an admission that such publications are prior art for patentability determination purposes.
Today's engineering systems sustain desirable performance by using well-designed control systems based on fundamental principles and mathematics. Many engineering breakthroughs and improvements in sensing and computation have helped to advance the field. Control systems currently play critical roles in many areas, including manufacturing, electronics, communications, transportation, computers, and networks, as well as many commercial and military systems. Traditionally, almost all modern control design is based on forcing the nonlinear systems to perform and behave like linear systems, thus limiting its maximum potential. In this paper a novel nonlinear control design methodology is introduced that overcomes this limitation.
Several of the popular advanced nonlinear control system approaches are based in passivity and dissipative control theories. Initially, P. J. Moylan, Implications of Passivity in a Class of Nonlinear Systems, IEEE Transactions of Automatic Control, Vol. AC-19, No. 4, August 1974, pp. 373-381, discussed the implications of passivity for a broad class of nonlinear systems, wherein a connection is established between the input-output property of passivity and a set of constraints on the state equations for the system. Later, J. L. Wyatt, Jr., L. O. Chua, J. W. Gannett, I. C. Goknar, and D. N. Green, Energy Concepts in the State-Space Theory of Nonlinear n-Ports: Part I—Passivity, IEEE Transactions on Circuits and Systems, Vol. CAS-28, No. 1, January 1981, pp. 48-61; and J. L. Wyatt, Jr., L. O. Chua, J. W. Gannett, I. C. Goknar, and D. N. Green, Energy Concepts in the State-Space Theory of Nonlinear n-Ports: Part II—Losslessness, IEEE Transactions on Circuits and Systems, Vol. CAS-29, No. 7, July 1982, pp. 417-430, clarified the meaning of passivity and losslessness as understood in nonlinear circuit theory, and their counterparts in classical physics. Most recently, R. Ortega, Z. P. Jiang, and D. J. Hill, Passivity-Based Control of Nonlinear systems: A Tutorial, Proceedings of the American Control Conference, Albuquerque, N. Mex., June 1997, pp. 2633-237, reviewed recent results on the stabilization of nonlinear systems using a passivity approach. Passivity properties play a vital role in designing asymptotically stabilizing controllers for nonlinear systems where the nonlinear versions of the Kalman-Yacubovitch-Popov lemma are used as key testing tools. The dissipative characteristics of dynamical systems have their origins in work by J. C. Willems, Dissipative Dynamical Systems Part I: General Theory; Part II: Linear Systems with Quadratic Supply Rates, Archive for Rational Mechanics and Analysis, Vol. 45, pp. 321-393, 1972, with further specifics given by D. Hill and P. J. Moylan, The Stability of Nonlinear Dissipative Systems, IEEE Transactions on Automatic Control, October 1976, pp. 708-710, in which a technique is introduced for generating Lyapunov functions for a broad class of nonlinear systems represented by state equations. The system, for which a Lyapunov function is required, is assumed to have a property called dissipativeness. In other words, the system absorbs more energy from the external world than it supplies. Different types of dissipativeness can be considered depending on how the “power input” is selected. Dissipativeness is shown to be characterized by the existence of a computable function which can be interpreted as the “stored energy” of the system. Under certain conditions, this energy function is a Lyapunov function which establishes stability, and in some cases asymptotic stability, of the isolated system. It was shown that for a certain class of nonlinear systems, that an “energy” approach was useful in analyzing stability. P. Kokotovic and M. Arcak, Constructive Nonlinear Control: A Historical Perspective, Preprint submitted to Elsevier, August 2000, provide a recent discussion about the historical perspective of constructive nonlinear control theories. Structural properties of nonlinear systems and passivation-based designs exploit the connections between passivity and inverse optimality, and between Lyapunov functions and optimal value functions. Recursive design procedures, such as backstepping and forwarding, achieve certain optimal properties for important classes of nonlinear systems. Some of the more popular nonlinear control system designs, e.g., J.-J. E. Slotine and W. Li, Applied Nonlinear Control, Prentice Hall, Inc., N.J., 1991; M. Kristic, I. Kanellakopoulos, and P. Kokotovic, Nonlinear and Adaptive Control Design, John Wiley & Sons, Inc., New York, 1995; and P. A. Ioannou and J. Sun, Robust Adaptive Control, Prentice Hall, Englewood Cliffs, N.J., 1996, have their fundamental foundations built upon these concepts.
In other engineering disciplines, A. A. Alonso and B. E. Ydstie, Stabilization of Distributed Systems Using Irreversible Thermodynamics, Automatica, Vol. 37, 2001, pp. 1739-1755, connect thermodynamics and the passivity theory of nonlinear control. The storage function is derived from the convexity of the entropy and is closely related to thermodynamic availability. Dissipation is related to positive entropy production. In this form the supply function is a product of force and flow variation variables. Results are discussed in relationship to heat conduction and reaction diffusion equation problems. K.-H. Anthony, Hamilton's Action Principle and Thermodynamics of Irreversible Processes—A Unifying Procedure for Reversible and Irreversible Processes, J. Non-Newtonian Fluid Mechanics, Vol. 96, 2001, pp. 291-339, suggests that non-equilibrium thermodynamics of irreversible processes may be included into the framework of a Lagrangian formalism. This formalism presents a unified method for reversible and irreversible processes. A straightforward procedure allows for the incorporation of both the first and second laws of thermodynamics into the Lagrangian. The theory is illustrated in three representative examples which include; material flow, heat conduction, diffusion and chemical reactions.
The main contribution of this paper is to present a novel nonlinear control design methodology that is based on thermodynamic exergy and irreversible entropy production concepts. Relationships are developed between exergy, irreversible entropy production, Hamiltonian systems, Lyapunov optimal functions, electric AC power concepts, and power flow, for control system design. Both necessary and sufficient conditions for stability are determined for nonlinear systems. By combining the first and second laws of thermodynamics, an exergy analysis approach is developed to construct Lyapunov optimal functions for Hamiltonian systems. The first time derivative of the Lyapunov functions, based on exergy, irreversible entropy production rate, and power flow is partitioned into either exergy dissipative or exergy generative terms.
Embodiments of the present invention are also described in Robinett, et al., “Exergy and Irreversible Entropy Production Thermodynamic Concepts for Control System Design: Robotic Servo Applications”, Proceedings of the 2006 IEEE International Conference on Robotics and Automation, pp. 3685-3692, May 2006; Robinett, et al., “Exergy and Irreversible Entropy Production Thermodynamic Concepts for Control System Design: Regulators”, Proceedings of the 2006 IEEE International Conference on Control Applications, pp. 2249-2256, October 2006; and Robinett, et al., “Exergy and Irreversible Entropy Production Thermodynamic Concepts for Control System Design: Nonlinear Systems”, Proceedings of the 2006 14th Mediterranean Conference on Control and Automation, pp. 1-8 (June 2006); Robinett III, R. D., Wilson, D. G. and Reed, A. W., “Exergy Sustainability for Complex Systems”, InterJournal Complex Systems, 1616, New England Complex Systems Institute, September 2006; Robinett III, R. D. and Wilson, D. G., “Decentralized Exergy/Entropy Thermodynamic Control for Collective Robotic Systems”, ASME 2007 International Mechanical Engineering Congress & Exposition, Seattle, Wash., Nov. 11-15, 2007; Robinett, III, R. D. and Wilson, D. G., “Collective Plume Tracing: A Minimal Information Approach to Collective Control,” IEEE 2007 American Control Conference, New York, N.Y., July 2007; Robinett III, R. D. and Wilson, D. G., “Exergy and Entropy Thermodynamic Concepts for Nonlinear Control Design”, 2006 ASME International Mechanical Engineering Congress & Exposition, Chicago, Ill., Nov. 5-10, 2006; Robinett III, R. D. and Wilson, D. G., “Exergy and Entropy Thermodynamic Concepts for Control System Design: Slewing Single Axis”, 2006 AIAA Guidance, Navigation, and Control Conference, Keystone, Colo., Aug. 21-24, 2006; Robinett III, R. D. and Wilson, D. G., “Exergy and Irreversible Entropy Production Thermodynamic Concepts for Control Design: Nonlinear Systems”, 14th Mediterranean Conference on Control and Automation, Ancona, Italy, Jun. 28-30, 2006; Robinett III, R. D. and Wilson, D. G., “Exergy and Irreversible Entropy Production Thermodynamic Concepts for Control Design: Nonlinear Regulator Systems”, The 8th IASTED International Conference on Control and Applications, Montreal, Quebec, Canada, May 24-26, 2006; Robinett III, R. D., Wilson, D. G., and Reed, A. W., “Exergy Sustainability”, Sandia Report SAND2006-2759, May 2006; and Robinett III, R. D. and Wilson, D. G., “Collective Systems: Physical and Information Exergies”, Sandia Report SAND2007-2327, April 2007, all of which are incorporated herein by reference.
BRIEF SUMMARY OF THE INVENTION
The present invention is of a control system design method and concomitant control system comprising: representing a physical apparatus to be controlled as a Hamiltonian system; determining elements of the Hamiltonian system representation which are power generators, power dissipators, and power storage devices; analyzing stability and performance of the Hamiltonian system based on the results of the determining step and determining necessary and sufficient conditions for stability of the Hamiltonian system; creating a stable control system based on the results of the analyzing step; and employing the resulting control system to control the physical apparatus. In the preferred embodiment, the control system is nonlinear. Stability is determined in terms of power flow, preferably in terms of whether the Hamiltonian system is moving toward or away from its minimum energy and maximum entropy state. Stability is determined via a Lyapunov derivative function, preferably one that is a decomposition and sum of exergy dissapation rate and exergy generation rate. The Hamiltonian system can be, for example, a forced harmonic oscillator, a mass-spring-damper dynamic system, or a Duffing oscillator/Coulomb friction nonlinear damper system. The physical apparatus can be any of the following: manufacturing devices, electronic devices, communications devices, transportation apparatuses, computer devices, computer networks, military apparatus, electrical power grids, Supervisory Control And Data Acquisition apparatus, and pipeline apparatus.
Further scope of applicability of the present invention will be set forth in part in the detailed description to follow, taken in conjunction with the accompanying drawings, and in part will become apparent to those skilled in the art upon examination of the following, or may be learned by practice of the invention. The objects and advantages of the invention may be realized and attained by means of the instrumentalities and combinations particularly pointed out in the appended claims.



BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS
The accompanying drawings, which are incorporated into and form a part of the specification, illustrate one or more embodiments of the present invention and, together with the description, serve to explain the principles of the invention. The drawings are only for the purpose of illustrating one or more preferred embodiments of the invention and are not to be construed as limiting the invention. In the drawings:
FIG. 1 illustrates the first law of thermodynamics for a volume control;
FIG. 2 illustrates entropy within a flux exchange system;
FIG. 3 shows time response for power in a general AC circuit with ω=2π, v=1.5, ī=2.0, and θ=π/4;
FIG. 4 illustrates a mass-spring-damper system;
FIG. 5 illustrates a neutrally stable reversible, conservative system;
FIGS. 6(a) and (b) show for case 1 entropy control harmonic oscillator results;
FIGS. 7(a) and (b) show for case 1a entropy control harmonic oscillator results;
FIGS. 8(a) and (b) show for case 2 entropy control harmonic oscillator results;
FIGS. 9(a) and (b) show for case 2a entropy control harmonic oscillator results;
FIGS. 10(a) and (b) show for case 3 entropy control harmonic oscillator results;
FIGS. 11(a) and (b) show for case 3a entropy control harmonic oscillator results;
FIG. 12 shows for case 1 Bode response at resonance results;
FIGS. 13(a) and (b) show for case 1 mass-spring-damper with PD control numerical results;
FIGS. 14(a) and (b) show for case 2 mass-spring-damper with PD control numerical results;
FIGS. 15(a) and (b) show for case 3 mass-spring-damper with PD control numerical results;
FIGS. 16(a) and (b) show for case 1a mass-spring-damper with PID control numerical results;
FIGS. 17(a) and (b) show for case 1 mass-spring-damper with PID control numerical results;
FIGS. 18(a) and (b) show for case 2 mass-spring-damper with PID control numerical results;
FIG. 19 shows for case 2a mass-spring-damper spring restoring and inertial effects numerical results;
FIGS. 20(a) and (b) show for case 3 mass-spring-damper with PID control numerical results;
FIGS. 21(a) and (b) show for all cases mass-spring-damper with PID control numerical results;
FIG. 22 illustrates a Duffing oscillator/Coulomb friction system;
FIGS. 23(a) and (b) show for case 1 Duffing oscillator/Coulomb friction with PID control numerical results;
FIGS. 24(a) and (b) show for case 2 Duffing oscillator/Coulomb friction with PID control numerical results;
FIGS. 25(a) and (b) show for case 3 Duffing oscillator/Coulomb friction with PID control numerical results;
FIGS. 26(a) and (b) show for case 4 Duffing oscillator/Coulomb friction with PID control numerical results;
FIGS. 27(a) and (b) show for all cases Duffing oscillator/Coulomb friction numerical results;
FIG. 28 illustrates gain scheduling with the integral gain as a function of initial conditions;
FIGS. 29(a) and (b) show for case 1 mass-spring-damper with PID tracking control numerical results;
FIGS. 30(a) and (b) show for case 1a mass-spring-damper with PID tracking control numerical results;
FIGS. 31(a) and (b) show for case 2 mass-spring-damper with PID tracking control numerical results;
FIGS. 32(a) and (b) show for case 2a mass-spring-damper with PID tracking control numerical results;
FIGS. 33(a) and (b) show for case 3 mass-spring-damper with PID tracking control numerical results; and
FIGS. 34(a) and (b) show for case 3a mass-spring-damper with PID tracking control numerical results.



DETAILED DESCRIPTION OF THE INVENTION
The present invention is of a method of designing control systems for physical apparatuses. For purposes of the specification and claims, physical apparatuses include devices used in the fields of robotics, automation, manufacturing, electronics, communications, transportation, computers, and computer networks.
The present invention is of a novel control system design methodology. The methodology provides a unique set of criteria to design nonlinear controllers for nonlinear systems with respect to both performance and stability. It works by combining: concepts from thermodynamic exergy and entropy; Hamiltonian systems; Lyapunov's direct method and Lyapunov optimal analysis; electric AC power concepts; and power flow analysis. The thermodynamic concepts are employed to allow the control design to be viewed as a power flow approach that is subject to the balance of power generation and power dissipation for the Hamiltonian system. This approach provides both necessary and sufficient conditions for stability while simultaneously allowing for performance specifications. Traditionally almost all modern control design is based on forcing the nonlinear systems to perform and behave like linear systems, thus limiting its maximum potential. In this invention a novel nonlinear control design methodology is introduced that overcomes this limitation.
The present invention can revolutionize the way that control systems are designed. The field of control plays a critical role in many areas, such as: manufacturing, electronics, communications, transportation, computers, and networks, as well as many commercial and military systems. This methodology is directly applicable to the design and control of complex, multi-component and often adaptive systems of arbitrary purpose, design and underlying physics, such as critical infrastructure (electric power grids, Supervisory Control And Data Acquisition (SCADA) systems, telecommunication and satellite systems, oil and gas pipelines) and military systems (network centric warfare, critical assets/nuclear security systems). The secure and reliable operation for all these systems is vital to a healthy economy for both civil and military infrastructure. In addition, the methodology is inclusive of exergy sustainability for most open systems which may include mass, energy, entropy, and exergy flow analysis. For example, this approach can be used to design efficient satellite spacecraft systems, distributed/decentralized energy system utilization, and vehicle survivability subject to maximizing finite power limitations during the course of the mission. This is a novel approach to the design of these types of systems.
The methodology is described and then is demonstrated with several fundamental numerical simulation examples: (1) a forced harmonic oscillator; (2) a mass-spring-damper dynamic system that employs PD/PID (proportional-derivative/proportional-integral-derivative) regulator control; (3) a Duffing oscillator/Coulomb friction nonlinear model that employs PID regulator control; and (4) a linear PID tracking control design for the mass-spring-damper system. The control system performances are partitioned and evaluated based on exergy generation and exergy dissipation terms. This novel nonlinear control methodology results in both necessary and sufficient conditions for stability of nonlinear systems.
In this section the first and second laws of thermodynamics are introduced. The first law of thermodynamics states energy is conserved (see FIG. 1). The second law of thermodynamics states that the entropy of the universe always increases. The first law is a conservation equation while the second law is an inequality. Mathematically, the first law can be written in terms of its time derivatives or energy rate for a system as




  
    
      
        
          E
          .
        
        =
        
          
            
              ∑
              i
              
                
              
            
            ⁢
            
              
            
            ⁢
            
              
                Q
                .
              
              i
            
          
          +
          
            
              ∑
              j
              
                
              
            
            ⁢
            
              
            
            ⁢
            
              
                W
                .
              
              j
            
          
          +
          
            
              ∑
              k
              
                
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  
                    m
                    .
                  
                  k
                
                ⁡
                
                  (
                  
                    
                      h
                      k
                    
                    +
                    
                      ke
                      k
                    
                    +
                    
                      pe
                      k
                    
                    +
                    ...
                  
                  )
                
              
              .
            
          
        
      
    
    
      
        (
        1
        )
      
    
  




The term on the left represents the rate at which energy is changing within the system. The heat entering or leaving the system is given by {dot over (Q)}i the work entering or leaving the system is given by {dot over (W)}j. Next, material can enter or leave the system by {dot over (m)}k that includes enthalpy, h, kinetic and potential energies, ke, pe, etc. In addition, each term is “summed” over an arbitrary number of entry and exit locations i, j, k.

The second law or entropy rate equation for a system is given as




  
    
      
        
          S
          .
        
        =
        
          
            
              
                ∑
                i
                
                  
                
              
              ⁢
              
                
                  
                    Q
                    .
                  
                  i
                
                
                  T
                  i
                
              
            
            +
            
              
                ∑
                k
                
                  
                
              
              ⁢
              
                
              
              ⁢
              
                
                  
                    m
                    .
                  
                  k
                
                ⁢
                
                  s
                  k
                
              
            
            +
            
              
                S
                .
              
              i
            
          
          =
          
            
              
                S
                .
              
              e
            
            +
            
              
                
                  S
                  .
                
                i
              
              .
            
          
        
      
    
    
      
        (
        2
        )
      
    
  




where the left hand term is the rate entropy changes within the system and the right hand terms represent, in order, the rate heat conducts entropy to and from the system and the rate material carries it in or out. These two terms can be combined into one term {dot over (S)}e, the entropy exchanged (either positive or negative) with the environment and {dot over (S)}i is the irreversible entropy production rate within the system. FIG. 2 shows the entropy exchanges and production within the system.

The irreversible entropy production rate can be written as the sum of the thermodynamic forces and the thermodynamic flows




  
    
      
        
          S
          .
        
        =
        
          
            
              ∑
              k
              
                
              
            
            ⁢
            
              
                F
                k
              
              ⁢
              
                
                  X
                  .
                
                k
              
            
          
          ≥
          0
        
      
    
    
      
        (
        3
        )
      
    
  




where the entropy change is the sum of all the changes due to the irreversible flows {dot over (X)}k with respect to each corresponding thermodynamic force Fk.

Next, for systems maintained at constant temperature and volume, a thermodynamic quantity called the Helmholtz free energy function is defined as

Ξ=E−ToS  (4)

where To is the reference environmental temperature. The Helmholtz free energy function is described as the maximum theoretically available energy that can do work which we call exergy. Exergy is also known as negative-entropy. By taking the time derivative of the Helmholtz free energy function (4) and substituting in the expressions for the first (1) and second (2) (multiplied by the environmental temperature) laws, results in the exergy rate equation





  
    
      
        
          Ξ
          .
        
        =
        
          
            
              ∑
              i
              
                
              
            
            ⁢
            
              
            
            ⁢
            
              
                (
                
                  1
                  -
                  
                    
                      T
                      o
                    
                    
                      T
                      i
                    
                  
                
                )
              
              ⁢
              
                
                  Q
                  .
                
                i
              
            
          
          +
          
            
              ∑
              j
              
                
              
            
            ⁢
            
              (
              
                
                  
                    W
                    .
                  
                  j
                
                -
                
                  
                    p
                    o
                  
                  ⁢
                  
                    
                      d
                      
                        V
                        _
                      
                    
                    
                      d
                      t
                    
                  
                
              
              )
            
          
          +
          
            
              ∑
              k
              
                
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  m
                  .
                
                k
              
              ⁢
              
                ζ
                k
                flow
              
            
          
          -
          
            
              T
              o
            
            ⁢
            
              
                
                  S
                  .
                
                i
              
              .
            
          
        
      
    
    
      
        (
        5
        )
      
    
  




where {dot over (Ξ)} is the rate at which exergy stored within the system is changing. The terms on the right, in order, define the rate exergy is carried in/out by; i) heat, ii) work (less any work the system does on the environment at constant environmental pressure po if the system volume V changes), and iii) by the material (or quantity known as flow exergy). The final term, To{dot over (S)}i, is the rate exergy is destroyed within the system.

Next, Hamiltonian Mechanics are discussed. The derivation of the Hamiltonian begins with the Lagrangian for a system defined as

L=T(q,{dot over (q)},t)−V(q,t)  (6)

where

    
    
        t=time explicitly
        q=N-dimensional generalized coordinate vector
        {dot over (q)}=N-dimensional generalized velocity vector
        T=Kinetic energy
        V=Potential energy
    
    


The Hamiltonian is defined in terms of the Lagrangian as




  
    
      
        
          H
          ≡
          
            
              
                ∑
                
                  i
                  =
                  1
                
                n
              
              ⁢
              
                
              
              ⁢
              
                
                  
                    ∂
                    L
                  
                  
                    ∂
                    
                      q
                      .
                    
                  
                
                ⁢
                
                  q
                  .
                
              
            
            -
            
              L
              ⁡
              
                (
                
                  q
                  ,
                  
                    q
                    .
                  
                  ,
                  t
                
                )
              
            
          
        
        =
        
          
            H
            ⁡
            
              (
              
                q
                ,
                
                  q
                  .
                
                ,
                t
              
              )
            
          
          .
        
      
    
    
      
        (
        7
        )
      
    
  




The Hamiltonian in terms of the canonical coordinates (q, p) is





  
    
      
        
          H
          ⁡
          
            (
            
              q
              ,
              p
              ,
              t
            
            )
          
        
        =
        
          
            
              ∑
              
                i
                =
                1
              
              n
            
            ⁢
            
              
            
            ⁢
            
              
                p
                i
              
              ⁢
              
                
                  q
                  .
                
                i
              
            
          
          -
          
            L
            ⁡
            
              (
              
                q
                ,
                
                  q
                  .
                
                ,
                t
              
              )
            
          
        
      
    
    
      
        (
        8
        )
      
    
  




where the canonical momentum is defined as





  
    
      
        
          p
          i
        
        =
        
          
            
              ∂
              L
            
            
              ∂
              
                q
                .
              
            
          
          .
        
      
    
    
      
        (
        9
        )
      
    
  




Then Hamilton's canonical equations of motion become





  
    
      
        
          
            
              
                
                  q
                  .
                
                i
              
              =
              
                
                  ∂
                  H
                
                
                  ∂
                  
                    p
                    i
                  
                
              
            
          
        
        
          
            
              
                
                  p
                  .
                
                i
              
              =
              
                
                  -
                  
                    
                      ∂
                      H
                    
                    
                      ∂
                      
                        q
                        i
                      
                    
                  
                
                +
                
                  Q
                  i
                
              
            
          
        
      
    
    
      
        (
        10
        )
      
    
  




where Qi is the generalized force vector. Next taking the time derivative of (8) gives





  
    
      
        
          H
          .
        
        =
        
          
            ∑
            
              i
              =
              1
            
            n
          
          ⁢
          
            
          
          ⁢
          
            
              (
              
                
                  
                    
                      p
                      .
                    
                    i
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    i
                  
                
                +
                
                  
                    p
                    i
                  
                  ⁢
                  
                    
                      q
                       ̈
                    
                    i
                  
                
                -
                
                  
                    ∂
                    L
                  
                  
                    ∂
                    t
                  
                
                -
                
                  
                    
                      ∂
                      L
                    
                    
                      ∂
                      
                        q
                        i
                      
                    
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    i
                  
                
                -
                
                  
                    
                      ∂
                      L
                    
                    
                      ∂
                      
                        
                          q
                          .
                        
                        i
                      
                    
                  
                  ⁢
                  
                    
                      q
                       ̈
                    
                    i
                  
                
              
              )
            
            .
          
        
      
    
    
      
        (
        11
        )
      
    
  




Then substituting (10) and simplifying gives





  
    
      
        
          H
          .
        
        =
        
          
            
              ∑
              
                i
                =
                1
              
              n
            
            ⁢
            
              
                Q
                i
              
              ⁢
              
                
                  q
                  .
                
                i
              
            
          
          -
          
            
              
                ∂
                L
              
              
                ∂
                t
              
            
            .
          
        
      
    
    
      
        (
        12
        )
      
    
  




Hamiltonians for most natural systems are not explicit functions of time (or ∂L/∂t=0). Then for

L=L(q,{dot over (q)})  (13)

the power (work/energy) equation becomes





  
    
      
        
          
            H
            .
          
          ⁡
          
            (
            
              q
              ,
              p
            
            )
          
        
        =
        
          
            ∑
            
              i
              =
              1
            
            n
          
          ⁢
          
            
              Q
              i
            
            ⁢
            
              
                
                  q
                  .
                
                i
              
              .
            
          
        
      
    
    
      
        (
        14
        )
      
    
  




Thermo-Mechanical Relationships are next discussed. A system is conservative if

{dot over (H)}=0 and H=constant.

A force is conservative if





  
        ⁢
    
      F
      ·
      dx
    
  
  =
  
    
            ⁢
      
        F
        ·
        vdt
      
    
    =
    
      
                ⁢
        
          Q
          j
        
        ⁢
        
          
            q
            .
          
          j
        
        ⁢
        dt
      
      =
      0
    
  




where F is the force, dx the displacement, and v the velocity. Basically, all of the forces can be modeled as potential force fields which are storage devices.

A thermodynamic system is reversible if




  
    
      
        dS
        =
        
        ⁢
        
          dQ
          T
        
      
    
  
  
    
      
        
                    ⁢
          
            d
            S
          
        
        =
        
        ⁢
        
          
                        ⁢
            
              
                d
                Q
              
              T
            
          
          =
          0
        
      
    
  
  
    
      
        
                    ⁢
          
            d
            S
          
        
        =
        
        ⁢
        
          
                        ⁡
            
              [
              
                
                  d
                  
                    S
                    i
                  
                
                +
                
                  d
                  
                    S
                    e
                  
                
              
              ]
            
          
          =
          
            
              
                                ⁡
                
                  [
                  
                    
                      
                        S
                        .
                      
                      i
                    
                    +
                    
                      
                        S
                        .
                      
                      e
                    
                  
                  ]
                
              
              ⁢
              
                d
                t
              
            
            =
            0
          
        
      
    
  




which implies that {dot over (S)}e={dot over (Q)}/T since by definition the second law gives {dot over (S)}i=0.

A thermodynamic system is irreversible if for

dS=[{dot over (S)}i+{dot over (S)}e]dt=0

then {dot over (S)}e≦0 and {dot over (S)}i≧0.

Now the connections between thermodynamics and Hamiltonian mechanics are investigated.
1. The irreversible entropy production rate can be expressed as




  
    
      
        
          
            S
            .
          
          i
        
        =
        
          
            
              ∑
              k
              
                
              
            
            ⁢
            
              
            
            ⁢
            
              
                F
                k
              
              ⁢
              
                
                  X
                  .
                
                k
              
            
          
          =
          
            
              
                1
                
                  T
                  o
                
              
              ⁢
              
                
                  ∑
                  k
                  
                    
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    Q
                    k
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    k
                  
                
              
            
            ≥
            0.
          
        
      
    
    
      
        (
        15
        )
      
    
  




2. The time derivative of the Hamiltonian is equivalent to the exergy rate




  
    
      
        
          
            
              
                H
                .
              
              =
              
              ⁢
              
                
                  ∑
                  k
                  
                    
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    Q
                    k
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    k
                  
                
              
            
          
        
        
          
            
              
                Ξ
                .
              
              =
              
              ⁢
              
                
                  
                    W
                    .
                  
                  -
                  
                    
                      T
                      o
                    
                    ⁢
                    
                      
                        S
                        .
                      
                      i
                    
                  
                
                =
                
                  
                    
                      ∑
                      
                        j
                        =
                        1
                      
                      N
                    
                    ⁢
                    
                      
                    
                    ⁢
                    
                      
                        Q
                        j
                      
                      ⁢
                      
                        
                          q
                          .
                        
                        j
                      
                    
                  
                  -
                  
                    
                      ∑
                      
                        l
                        =
                        N
                      
                      
                        M
                        -
                        N
                      
                    
                    ⁢
                    
                      
                    
                    ⁢
                    
                      
                        Q
                        l
                      
                      ⁢
                      
                        
                          
                            q
                            .
                          
                          l
                        
                        .
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        16
        )
      
    
  




Where the following assumptions apply when utilizing the exergy rate equation (5) for Hamiltonian systems:
(a) No substantial heat flow:

{dot over (Q)}i≈0.

(b) No substantial exergy flow or assume Ti is only slightly greater than To:




  
    1
    -
    
      
        T
        o
      
      
        T
        i
      
    
  
  ≈
  0.




(c) No po V work on the environment:




  
    
      p
      o
    
    ⁢
    
      
        d
        
          V
          _
        
      
      
        d
        t
      
    
  
  =
  0.




(d) No mass flow rate:




  
    
      ∑
      k
      
        
      
    
    ⁢
    
      
    
    ⁢
    
      
        
          m
          .
        
        k
      
      ⁢
      
        ζ
        k
        flow
      
    
  
  =
  0.




(e) Then define:

{dot over (W)}≧0 powerinput/generated

To{dot over (S)}i≧0 powerdissipated.

3. A conservative system is equivalent to a reversible system when

{dot over (H)}=0 and {dot over (S)}e=0

then

{dot over (S)}i=0 and {dot over (W)}=0.

4. For a system that “appears to be conservative”, but is not reversible is defined as:




  
    
      
        
          
            H
            .
          
          ave
        
        =
        
        ⁢
        0
      
    
    
      
        =
        
        ⁢
        
          
            1
            τ
          
          ⁢
          
                        ⁡
            
              [
              
                
                  W
                  .
                
                -
                
                  
                    T
                    o
                  
                  ⁢
                  
                    
                      S
                      .
                    
                    i
                  
                
              
              ]
            
          
          ⁢
          
            d
            t
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            
              (
              
                W
                .
              
              )
            
            ave
          
          -
          
            
              (
              
                
                  T
                  o
                
                ⁢
                
                  
                    S
                    .
                  
                  i
                
              
              )
            
            ave
          
        
      
    
    
      
        =
        
        ⁢
        
          
            1
            τ
          
          ⁢
          
                        ⁡
            
              [
              
                
                  
                    ∑
                    j
                    N
                  
                  ⁢
                  
                    
                  
                  ⁢
                  
                    
                      Q
                      j
                    
                    ⁢
                    
                      
                        q
                        .
                      
                      j
                    
                  
                
                -
                
                  
                    ∑
                    
                      l
                      =
                      N
                    
                    
                      M
                      -
                      N
                    
                  
                  ⁢
                  
                    
                  
                  ⁢
                  
                    
                      Q
                      l
                    
                    ⁢
                    
                      
                        q
                        .
                      
                      l
                    
                  
                
              
              ]
            
          
          ⁢
          
            d
            t
          
        
      
    
  
  
    
      
        =
        
        ⁢
        averagepoweroveracycle
      
    
    
      
        
      
    
  




where τ is the period of the cycle. To be more specific about the average power calculations, the AC power factor provides an excellent example. For the general case of alternating current supplied to a complex impedance the voltage and current differ in phase by an angle θ. The time responses for power, voltage, and current are shown for a general AC circuit in FIG. 3. For





  
    
      
        
          W
          .
        
        =
        
        ⁢
        
          P
          =
          
            
              Q
              ⁢
              
                q
                .
              
            
            =
            
              vi
              =
              
                
                  2
                
                ⁢
                
                  v
                  _
                
                ⁢
                
                  
                    cos
                    ⁡
                    
                      (
                      
                        
                          ω
                          ⁢
                          
                            
                          
                          ⁢
                          t
                        
                        +
                        θ
                      
                      )
                    
                  
                  ·
                  
                    2
                  
                
                ⁢
                
                  i
                  _
                
                ⁢
                cos
                ⁢
                
                  
                
                ⁢
                ω
                ⁢
                
                  
                
                ⁢
                t
              
            
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            v
            _
          
          ⁢
          
            
              i
              _
            
            ⁡
            
              [
              
                
                  cos
                  ⁢
                  
                    
                  
                  ⁢
                  θ
                
                +
                
                  cos
                  ⁡
                  
                    (
                    
                      
                        2
                        ⁢
                        ω
                        ⁢
                        
                          
                        
                        ⁢
                        t
                      
                      +
                      θ
                    
                    )
                  
                
              
              ]
            
          
        
      
    
  




where P is power, v is voltage ( v), i is current (ī), θ is the phase angle, and ω is the frequency. Integrating over a cycle gives

({dot over (W)})ave= vī cos θ

where for the second term

 cos(2ωt+θ)dt=0.

This is an important set of conditions that will be used in the next section to find the generalized stability boundary.
5. Finally, the power terms are sorted into three categories:
(a) ({dot over (W)})ave−power generators; (Qj{dot over (q)}j)ave>0
(b) (To{dot over (S)}i)ave−power dissipators; (Qi{dot over (q)}i)ave<0
(c) (Qk {dot over (q)}k)ave=0; reversible/conservative exergy storage terms.
These three categories are fundamental terms in the following design procedures.
Next discussed are the necessary and sufficient conditions for stability. The Lyapunov function is defined as the total energy which for most mechanical systems is equivalent to an appropriate Hamiltonian function

V=H  (17)

which is positive definite. The time derivative is





  
    
      
        
          
            
              
                V
                .
              
              =
              
              ⁢
              
                
                  H
                  .
                
                =
                
                  
                    
                      ∑
                      k
                      
                        
                      
                    
                    ⁢
                    
                      
                    
                    ⁢
                    
                      
                        Q
                        k
                      
                      ⁢
                      
                        
                          q
                          .
                        
                        k
                      
                    
                  
                  =
                  
                    
                      
                        ∑
                        j
                        N
                      
                      ⁢
                      
                        
                      
                      ⁢
                      
                        
                          Q
                          j
                        
                        ⁢
                        
                          
                            q
                            .
                          
                          j
                        
                      
                    
                    -
                    
                      
                        ∑
                        
                          l
                          =
                          N
                        
                        
                          M
                          -
                          N
                        
                      
                      ⁢
                      
                        
                      
                      ⁢
                      
                        
                          Q
                          l
                        
                        ⁢
                        
                          
                            q
                            .
                          
                          l
                        
                      
                    
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                
                  W
                  .
                
                -
                
                  
                    T
                    o
                  
                  ⁢
                  
                    
                      
                        S
                        .
                      
                      i
                    
                    .
                  
                
              
            
          
        
      
    
    
      
        (
        18
        )
      
    
  




To describe a nonlinear system's behavior, two theorems help to characterize the essential features of its motion. In addition, by bounding the Lyapunov function between these Theorems, both necessary and sufficient conditions are a result of the transition of the time derivative of the Lyapunov function from stable to unstable.

    
    
        1. Lyapunov Theorem for Stability. Assume that there exists a scalar function V of the state x, with continuous first order derivatives such that

V(x) ispositivedefinite

{dot over (V)}(x) isnegativedefinite

V(x)→∞ as ∥x∥→∞

    
    


Then the equilibrium at the origin is globally asymptotically stable.

    
    
        2. Chetaev Theorem for Instability. Considering the equations of disturbed motion, let V be zero on the boundary of a region R which has the origin as a boundary point, and let both V and {dot over (V)} be positive-definite in R; then the undisturbed motion is unstable at the origin.
    
    


Based on the relationship between thermodynamic exergy and Hamiltonian systems a fundamental stability Lemma can be formulated: The stability of Hamiltonian systems is bounded between Theorems 1 and 2. Given the Lyapunov derivative as a decomposition and sum of exergy dissipation rate and exergy generation rate then:




  
    
      
        
          V
          .
        
        =
        
          
            
              W
              .
            
            -
            
              
                T
                o
              
              ⁢
              
                
                  S
                  .
                
                i
              
            
          
          =
          
            
              
                ∑
                
                  j
                  =
                  1
                
                N
              
              ⁢
              
                
              
              ⁢
              
                
                  Q
                  j
                
                ⁢
                
                  
                    q
                    .
                  
                  j
                
              
            
            -
            
              
                ∑
                
                  l
                  =
                  N
                
                
                  M
                  -
                  N
                
              
              ⁢
              
                
              
              ⁢
              
                
                  Q
                  l
                
                ⁢
                
                  
                    q
                    .
                  
                  l
                
              
            
          
        
      
    
    
      
        (
        19
        )
      
    
  




that is subject to the following general necessary and sufficient conditions:

To{dot over (S)}i≧0 Positivesemi-definite, alwaystrue

{dot over (W)}≧0 Positivesemi-definite; exergypumpedintothesystem.

The following corollaries encompass both stability and instability for Hamiltonian systems which utilize AC power concepts:

    
    
        Corollary 1: For (To{dot over (S)}i)ave=0 and ({dot over (W)})ave=0 then {dot over (V)}=0 the Hamiltonian system is neutrally stable, conservative and reversible.
        Corollary 2: For (To{dot over (S)}i)ave=0 and ({dot over (W)})ave>0 then {dot over (V)}>0 the Hamiltonian system is unstable.
        Corollary 3: For (To{dot over (S)}i)ave>0 and ({dot over (W)})ave=0 then {dot over (V)}<0 the Hamiltonian system is asymptotically stable and a passive system in the general sense (passivity controllers).
        Corollary 4: Given apriori (To{dot over (S)}i)ave>0 and ({dot over (W)})ave>0 then the Hamiltonian system is further subdivided into:
        
            4.1: For (To{dot over (S)}i)ave>({dot over (W)})ave with {dot over (V)}<0 yields asymptotic stability
            4.2: For (To{dot over (S)}i)ave=({dot over (W)})ave with {dot over (V)}=0 yields neutral stability
            4.3: For (To{dot over (S)}i)ave<({dot over (W)})ave with {dot over (V)}>0 yields an unstable system.

The bottom line is that stability is defined in terms of power flow which determines whether the system is moving toward or away from its minimum energy and maximum entropy state.

        
        
    
    


Present day applications use feedback controller designs that are Lyapunov Optimal. A control law is Lyapunov Optimal if it minimizes the first time derivative of the Lyapunov function over a space of admissible controls. In general, a set of feedback gains are optimized by minimizing the regulating and/or tracking error of the feedback controller while regulating to zero and/or tracking a desired reference input. The Lyapunov function is the total error energy which for most mechanical systems is equivalent to an appropriate Hamiltonian function

V=H.  (20)

Then the concept of Lyapunov Optimal follows directly from setting {dot over (W)}=0 in (19) and maximizing To{dot over (S)}i for which the time derivative of the Lyapunov function (Hamiltonian) or the modified power (work/energy) equation is written as





  
    
      
        
          V
          .
        
        =
        
          
            H
            .
          
          =
          
            
              
                -
                
                  T
                  o
                
              
              ⁢
              
                
                  S
                  .
                
                i
              
            
            =
            
              
                -
                
                  
                    ∑
                    
                      j
                      =
                      1
                    
                    N
                  
                  ⁢
                  
                    
                  
                  ⁢
                  
                    
                      Q
                      j
                    
                    ⁢
                    
                      
                        q
                        .
                      
                      j
                    
                  
                
              
              =
              
                -
                
                  
                    ∑
                    
                      j
                      =
                      1
                    
                    N
                  
                  ⁢
                  
                    
                  
                  ⁢
                  
                    
                      F
                      j
                    
                    ⁢
                    
                      
                        R
                        .
                      
                      j
                    
                  
                
              
            
          
        
      
    
    
      
        (
        21
        )
      
    
  




which is independent of system dynamics and is a kinematic quantity that applies to any system. Note that Fj denotes a set of forces acting on a mechanical system and {dot over (R)}j denotes the inertial linear velocity of the point where Fj is applied.

Passivity control follows directly from setting {dot over (W)}=0 in (19).
To summarize, the present invention provides a novel control system design methodology that uniquely combines: concepts from thermodynamic exergy and entropy; Hamiltonian systems; Lyapunov's direct method and Lyapunov optimal analysis; electric AC power concepts; and power flow analysis. Relationships are derived between exergy/entropy and Lyapunov optimal functions for Hamiltonian systems. This methodology is demonstrated below with several fundamental numerical simulation examples: 1) a forced harmonic oscillator; 2) a mass-spring-damper dynamic system that employs PD/PID regulator control; 3) a Duffing oscillator/Coulomb friction nonlinear model that employs PID regulator control; and 4) a linear PID tracking control design for the mass-spring-damper system. The control system performance results were partitioned and evaluated based on exergy generation and exergy dissipation terms. These results showed the stability boundaries for each numerical example which included both linear and nonlinear systems. This novel nonlinear control methodology resulted in both necessary and sufficient conditions for stability of nonlinear systems. In the near future, this novel control system design methodology will be extended to tracking and adaptive control of nonlinear systems. Although the methodology was demonstrated on simple linear and nonlinear dynamic systems, the methodology is applicable to a large class of nonlinear systems.
INDUSTRIAL APPLICABILITY
The invention is further illustrated by the following non-limiting examples. Several simple dynamic systems are investigated to demonstrate the exergy/entropy control design method of the invention.
EXAMPLE 1
The dynamic equation of motion for a forced harmonic oscillator (see FIG. 4) is defined as

M{umlaut over (x)}+C{umlaut over (x)}+Kx=Fo cos Ωt.  (22)

The steady-state response has the form

xp=X cos(Ωt−φ)  (23)

Substitute (23) into (22) and solve for X which produces





  
    X
    =
    
      
        
          
            
              F
              o
            
            /
            K
          
          
            
              
                
                  [
                  
                    1
                    -
                    
                      
                        (
                        
                          Ω
                          
                            ω
                            n
                          
                        
                        )
                      
                      2
                    
                  
                  ]
                
                2
              
              +
              
                
                  [
                  
                    2
                    ⁢
                    
                      ζ
                      ⁡
                      
                        (
                        
                          Ω
                          
                            ω
                            n
                          
                        
                        )
                      
                    
                  
                  ]
                
                2
              
            
          
        
        ⁢
        
          
        
        ⁢
        where
        ⁢
        
          
        
        ⁢
        
          ω
          n
          2
        
      
      =
      
        K
        M
      
    
  
  ,
  
    
  
  ⁢
  
    ζ
    =
    
      
        
          C
          /
          2
        
        ⁢
        M
        ⁢
        
          
        
        ⁢
        
          ω
          n
        
        ⁢
        
          
        
        ⁢
        and
        ⁢
        
          
        
        ⁢
        tan
        ⁢
        
          
        
        ⁢
        φ
      
      =
      
        2
        ⁢
        ζ
        ⁢
        
          
        
        ⁢
        
          (
          
            Ω
            
              ω
              n
            
          
          )
        
        ⁢
        
          
            1
            
              [
              
                1
                -
                
                  
                    (
                    
                      Ω
                      
                        ω
                        n
                      
                    
                    )
                  
                  2
                
              
              ]
            
          
          .
        
      
    
  




Also note that for a static deflection Xo=Fo/K, where the dynamic magnification factor becomes





  
    H
    ⁡
    
      (
      Ω
      )
    
  
  =
  
    
      X
      
        X
        o
      
    
    =
    
      
        
          1
          
            1
            -
            
              r
              2
            
          
        
        ⁢
        
          
        
        ⁢
        where
        ⁢
        
          
        
        ⁢
        r
      
      =
      
        
          
            Ω
            
              ω
              n
            
          
          ⁢
          
            
          
          ⁢
          for
          ⁢
          
            
          
          ⁢
          ζ
        
        =
        0.
      
    
  




For r=1 or Ω=ωn, the response is undefined or unstable. This is shown conventionally with the Bode plot (see FIG. 12).

The forced harmonic oscillator dynamic model instability can also be predicted with exergy/entropy control. Applying the exergy/entropy control design (for power input or exergy generation rate only) gives




  
    
      
        
          V
          .
        
        =
        
        ⁢
        
          
            W
            .
          
          -
          
            
              T
              o
            
            ⁢
            
              
                S
                .
              
              i
            
          
        
      
    
    
      
        =
        
        ⁢
        
          
            
              ∑
              
                j
                =
                1
              
              N
            
            ⁢
            
              
            
            ⁢
            
              
                Q
                j
              
              ⁢
              
                
                  q
                  .
                
                j
              
            
          
          -
          
            
              ∑
              
                l
                =
                N
              
              
                M
                -
                N
              
            
            ⁢
            
              
            
            ⁢
            
              
                Q
                l
              
              ⁢
              
                
                  q
                  .
                
                l
              
            
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          W
          .
        
      
    
    
      
        =
        
        ⁢
        
          Q
          ⁢
          
            
              q
              .
            
            .
          
        
      
    
  




Then by applying the following definitions for Q and {dot over (q)} at Ω=ωn, φ=π/2, and ζ=0 gives

Q=Fo cos Ωt 

{dot over (q)}=−XΩ sin(Ωt−π/2)

and with the appropriate variable substitutions yields





  
    
      
        
          Q
          ⁢
          
            q
            .
          
        
        =
        
        ⁢
        
          
            -
            
              F
              o
            
          
          ⁢
          X
          ⁢
          
            
          
          ⁢
          ΩcosΩ
          ⁢
          
            
          
          ⁢
          t
          ⁢
          
            
          
          ⁢
          
            sin
            ⁡
            
              (
              
                
                  Ω
                  ⁢
                  
                    
                  
                  ⁢
                  t
                
                -
                
                  π
                  /
                  2
                
              
              )
            
          
        
      
    
  
  
    
      
        
          =
          
          ⁢
          
            
              1
              2
            
            ⁢
            
              F
              o
            
            ⁢
            X
            ⁢
            
              
            
            ⁢
            
              Ω
              ⁡
              
                [
                
                  
                    sin
                    ⁡
                    
                      (
                      
                        π
                        /
                        2
                      
                      )
                    
                  
                  -
                  
                    sin
                    ⁡
                    
                      (
                      
                        
                          2
                          ⁢
                          Ω
                          ⁢
                          
                            
                          
                          ⁢
                          t
                        
                        -
                        
                          π
                          /
                          2
                        
                      
                      )
                    
                  
                
                ]
              
            
            ⁢
            
              
            
            ⁢
            where
          
        
        ⁢
        
          

        
        ⁢
        
          
        
        ⁢
        
          
            
              (
              
                W
                .
              
              )
            
            ave
          
          =
          
            
              
                (
                
                  Q
                  ⁢
                  
                    q
                    .
                  
                
                )
              
              ave
            
            =
            
              
                
                  1
                  2
                
                ⁢
                
                  F
                  o
                
                ⁢
                X
                ⁢
                
                  
                
                ⁢
                Ω
              
              >
              0.
            
          
        
      
    
  




The system is frequency and phase-locked which shows unbounded exergy generation and the response becomes infinite and unstable. This is also predicted in the following numerical simulations (see Case 2) with exergy/entropy control where the response grows unbounded.

Given the SDOF harmonic oscillator with zero inputs as

M{umlaut over (x)}+Kx=0

that is subject to the initial condition

x(0)=xo=1.0

then one can apply the entropy control design as





  
    
      
        
          V
          .
        
        =
        
        ⁢
        
          
            W
            .
          
          -
          
            
              T
              o
            
            ⁢
            
              
                S
                .
              
              i
            
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            
              ∑
              
                j
                =
                1
              
              N
            
            ⁢
            
              
            
            ⁢
            
              
                Q
                j
              
              ⁢
              
                
                  q
                  .
                
                j
              
            
          
          -
          
            
              ∑
              
                l
                =
                N
              
              
                M
                -
                N
              
            
            ⁢
            
              
            
            ⁢
            
              
                Q
                l
              
              ⁢
              
                
                  q
                  .
                
                l
              
            
          
        
      
    
  




which yields

(To{dot over (S)}i)ave=0 ({dot over (W)})ave=0 (To{dot over (S)}rev)ave=([M{umlaut over (x)}·{dot over (x)}+Kx·{dot over (x)}])ave=0.

This demonstrates Corollary 1 for a neutrally stable, reversible, conservative Hamiltonian system. Note that the forces in the system are conservative, such that the potential and kinetic energy rates or average powers are zero over each cycle (see FIG. 5).

Next, equation (22) with C=0 is explored with an excitation force. The corresponding time derivative of Lyapunov function/Hamiltonian is partitioned as a sum of exergy generation and exergy dissipation as




  
    V
    .
  
  =
  
    
      H
      .
    
    =
    
      
        
          W
          .
        
        -
        
          
            T
            o
          
          ⁢
          
            
              S
              .
            
            i
          
        
      
      =
      
        
          
            ∑
            
              J
              =
              1
            
            N
          
          ⁢
          
            
          
          ⁢
          
            
              Q
              j
            
            ⁢
            
              
                q
                .
              
              j
            
          
        
        -
        
          
            ∑
            
              l
              =
              N
            
            
              M
              -
              N
            
          
          ⁢
          
            
          
          ⁢
          
            
              Q
              l
            
            ⁢
            
              
                q
                .
              
              l
            
          
        
      
    
  




which yields

To{dot over (S)}i=0

{dot over (W)}=Fo cos Ωt·{dot over (x)}

(To{dot over (S)}rev)ave=(M{umlaut over (x)}·{dot over (x)}+Kx·{dot over (x)})ave=0.

Numerical simulations are performed for three separate cases with the numerical values listed in Table 1. These include before, at, and after resonance. The harmonic oscillator is subject to the same initial condition as previously. FIGS. 6(a) and (b) plot the input force with the velocity and the exergy and exergy rate for case 1, respectively. FIGS. 7(a) and (b) plot the system responses for position, velocity, and acceleration and the potential and kinetic energy rates for case 1, respectively. The beat frequency cycle is defined at approximately 25 seconds as can be observed by the zero crossings. For case 2, at resonance, the same responses are given and shown in FIGS. 8(a) and (b) and 9(a) and (b), respectively. Note that for a constant amplitude sinusoidal input force the velocity resonance grows asymptotically and unstable along with the system, potential/kinetic energy rates, and exergy/exergy rate responses. The exergy response grows accumulatively. For case 3 corresponding plots are given in FIGS. 10(a) and (b) and 11(a) and (b), respectively. For this case it is after resonance and the system oscillates with an accumulative affect being nulled out at each natural frequency cycle zero crossing, for potential/kinetic rates and exergy/exergy rate responses.







TABLE 1







Mass-spring harmonic oscillator cases











Case
Ω
Fo,
M
K


No.
(rad/sec)
(N)
(kg) 
(N/m)





1
 0.75
10.0
10.0
10.0


2
1.0
10.0
10.0 
10.0


3
5.0
10.0
10.0 
10.0









At resonance case 2 demonstrates Corollary 2 such that

(To{dot over (S)}i)ave=0 and ({dot over (W)})ave>0.  (24)

Case 2 also numerically demonstrated the instability of the harmonic oscillator at resonance which corresponds to the previous analytical analysis and conventional Bode plot (see FIG. 12). Cases one and three are before and after resonance and demonstrate over the proper cycle, the RMS values are ({dot over (W)})ave=0. This demonstrates Corollary 1 and reproduces the results of a linear system analyses.

EXAMPLE 2
The mass-spring-damper dynamic model (see FIG. 4) is defined as

M{umlaut over (x)}+C{dot over (x)}+Kx=u  (25)

where M, C, K, and u are the mass, damper, stiffness coefficients and external force input terms. The Proportional-Derivative (PD) and Proportional-Integral-Derivative (PID) controller is defined as





  
    
      
        u
        =
        
          
            
              -
              
                K
                P
              
            
            ⁢
            x
          
          -
          
            
              K
              I
            
            ⁢
            
              
                ∫
                0
                t
              
              ⁢
              
                x
                ⁢
                
                  
                
                ⁢
                
                  d
                  τ
                
              
            
          
          -
          
            
              K
              D
            
            ⁢
            
              x
              .
            
          
        
      
    
    
      
        (
        26
        )
      
    
  




where Kp, KI, and KD are the proportional, integral and derivative controller gains. Also note that for PD control, KI=0.0.

The passive PD control law is partitioned into terms of exergy dissipation and exergy generation. Now apply exergy/entropy control design and the time derivative of the Lyapunov function/Hamiltonian becomes




  
    
      
        
          V
          .
        
        =
        
          
            H
            .
          
          =
          
            
              
                W
                .
              
              -
              
                
                  T
                  o
                
                ⁢
                
                  
                    S
                    .
                  
                  i
                
              
            
            =
            
              
                
                  ∑
                  
                    j
                    =
                    1
                  
                  N
                
                ⁢
                
                  
                
                ⁢
                
                  
                    Q
                    j
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    j
                  
                
              
              -
              
                
                  ∑
                  
                    l
                    =
                    N
                  
                  
                    M
                    -
                    N
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    Q
                    l
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    l
                  
                
              
            
          
        
      
    
    
      
        (
        27
        )
      
    
  




which yields

To{dot over (S)}i=(C+KD){dot over (x)}·{dot over (x)}

{dot over (W)}=0

(To{dot over (S)}rev)ave=(M{umlaut over (x)}·{dot over (x)}+(K+KP)x·{dot over (x)})ave=0.  (28)

The first term in (28) is identified as a dissipative term composed of the derivative control term along with the damping term. The second term in (28) is zero or no generative terms present, characteristic of any PD control system. The proportional control term is added to the system stiffness term and contributes to the reversible portion or third term in (28). Hence the combined dynamic and PD control system is categorized as a passively stable system.
Numerical simulations are performed for three separate PD control regulator cases with the numerical values listed in Table 2. These include three different levels of damping. The mass-spring-damper system is subject to an initial condition of x0=1.0. The position, velocity, and acceleration responses along with the exergy and exergy rate responses are given in FIGS. 13(a) and (b), 14(a) and (b), and 15(a) and (b) for cases 1, 2, and 3, respectively. In all cases, the exergy rate is negative or dissipative and all the responses are stable and decaying as one would expect for PD passivity control.







TABLE 2







Mass-spring-damper model and PD/PID control system gains













Case
KP
KI(PD)/KI(PID)
KD(PD)/KD(PID)
M
K
C


No.
(kg/s2)
(kg/s3)
(kg/s)
(kg)
(kg/s2)
(kg/s)
















1A
10.0
NA/1.0
NA/0.0
10.0
10.0
0.0


1
10.0
0.0/1.0
1.6/2.0
10.0
10.0
0.1


2
10.0
0.0/1.0
0.0/0.4
10.0
10.0
0.1


3
10.0
0.0/1.0
6.4/0.0
10.0
10.0
0.1









All three cases demonstrate Corollary 3 such that

(To{dot over (S)}i)ave>0 and {dot over (W)}=0.  (29)

It is well known that PD control has been proven to be stable based on passivity control theory. Exergy/entropy control design has shown the PD controller to be passive stable and that conventional passivity control theory is a subset of this general development.

The PID control law is partitioned into terms of exergy generation and exergy dissipation. Now apply exergy/entropy control design and the time derivative of the Lyapunov function/Hamiltonian becomes




  
    
      
        
          V
          .
        
        =
        
          
            H
            .
          
          =
          
            
              
                W
                .
              
              -
              
                
                  T
                  o
                
                ⁢
                
                  
                    S
                    .
                  
                  i
                
              
            
            =
            
              
                
                  ∑
                  
                    j
                    =
                    1
                  
                  N
                
                ⁢
                
                  
                
                ⁢
                
                  
                    Q
                    j
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    j
                  
                
              
              -
              
                
                  ∑
                  
                    l
                    =
                    N
                  
                  
                    M
                    -
                    N
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    Q
                    l
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    l
                  
                
              
            
          
        
      
    
    
      
        (
        30
        )
      
    
  




which yields





  
    
      
        
          
            
              
                
                  T
                  o
                
                ⁢
                
                  
                    S
                    .
                  
                  i
                
              
              =
              
              ⁢
              
                
                  (
                  
                    C
                    +
                    
                      K
                      D
                    
                  
                  )
                
                ⁢
                
                  
                    x
                    .
                  
                  ·
                  
                    x
                    .
                  
                
              
            
          
        
        
          
            
              
                W
                .
              
              =
              
              ⁢
              
                
                  -
                  
                    K
                    I
                  
                
                ⁢
                
                  
                    ∫
                    0
                    t
                  
                  ⁢
                  
                    x
                    ⁢
                    
                      
                    
                    ⁢
                    
                      
                        d
                        τ
                      
                      ·
                      
                        x
                        .
                      
                    
                  
                
              
            
          
        
        
          
            
              
                
                  (
                  
                    
                      T
                      o
                    
                    ⁢
                    
                      
                        S
                        .
                      
                      rev
                    
                  
                  )
                
                ave
              
              =
              
              ⁢
              
                
                  
                    (
                    
                      
                        M
                        ⁢
                        
                          
                            x
                             ̈
                          
                          ·
                          
                            x
                            .
                          
                        
                      
                      +
                      
                        
                          (
                          
                            K
                            +
                            
                              K
                              P
                            
                          
                          )
                        
                        ⁢
                        
                          x
                          ·
                          
                            x
                            .
                          
                        
                      
                    
                    )
                  
                  ave
                
                =
                0.
              
            
          
        
      
    
    
      
        (
        31
        )
      
    
  




Once again, the first term in (31) is identified as a dissipative term composed of the derivative control term along with the damping term. The second term in (31) is the integral action and is identified as the generative input. Note that the proportional control term is added to the system stiffness term and contributes to the reversible portion of the third term in (31).

For PID control the stability boundary is determined as

{dot over (V)}=0={dot over (W)}−To{dot over (S)}i ∀t0≦t≦tfinal 

and is also considered a special case of the average power. In general for nonlinear systems, the average power terms will need to be taken into account.

Numerical simulations are performed for four separate PID control regulator cases with the numerical values listed in Table 2. The mass-spring-damper system is subject to an initial condition of x0=1.0. For case 1A, pure generative input is demonstrated where growing state oscillations (see FIG. 16(a)) and positively increasing values are observed in the responses for exergy and exergy rate (see FIG. 16(b)). The system has experienced pure input without any dissipation mechanisms to damp out the exergy coming into the system and looks very similar to the forced harmonic oscillator at resonance. For case 1 the integral of position, position, velocity, and acceleration responses along with the exergy and exergy rate responses are plotted in FIGS. 17(a) and (b). For this case, the dissipative term is greater than the generative term. This is observed from the decaying system responses. In case 2 the system responses along with the exergy and exergy rate responses are shown in FIGS. 18(a) and (b). In this case, the dissipative term is equal to the generative term. This results in system responses that do not decay or have constant oscillatory behavior. Further investigation of the exergy and exergy rate responses show that the mirror images for dissipation and generation will cancel each other out. For case 2 the potential and kinetic energy rate responses (see FIGS. 19(a) and (b)) are shown to have zero contribution or are reversible over each cycle. In case 3 the system responses along with the exergy and exergy rate responses are shown in FIGS. 20(a) and (b). In this case, the dissipative term is less than the generative term. This results in system responses with increasing oscillatory behavior. In addition, the sum of generative and dissipative terms are increasingly positive, hence exergy is being pumped into the system at a greater rate than can be dissipated. In conclusion, FIGS. 21(a) and (b) show the responses for the total exergy and exergy rates with respect to each case. Again, for case 2 a balance or boundary (neutral stability) can be observed.
The three cases for the PID control regulator mass-spring-damper dynamic system demonstrate the three subcases for Corollary 4: Given apriori To{dot over (S)}i>0 and {dot over (W)}>0 then the mass-spring-damper with PID control showed the following:

    
    
        4.1 Case 1 yielded (To{dot over (S)}i)ave>({dot over (W)})ave; asymptotic stability; damped stable response; and demonstration of Corollary 4.1.
        4.2 Case 2 yielded (To{dot over (S)}i)ave=({dot over (W)})ave; neutral stability; and demonstration of Corollary 4.2. This case is the dividing line where derivative and integral action cancel each other out.
        4.3 Case 3 yielded (To{dot over (S)}i)ave<({dot over (W)})ave; increasing unstable system; and demonstration of Corollary 4.3.
    
    


For PID control validation, first, compare a PID controller




  
    
      
        
          
            M
            ⁢
            
              x
               ̈
            
          
          +
          
            C
            ⁢
            
              x
              .
            
          
          +
          Kx
        
        =
        
          
            u
            PID
          
          =
          
            
              
                -
                
                  K
                  P
                
              
              ⁢
              x
            
            -
            
              
                K
                I
              
              ⁢
              
                
                  ∫
                  o
                  t
                
                ⁢
                
                  x
                  ⁢
                  
                    
                  
                  ⁢
                  
                    d
                    τ
                  
                
              
            
            -
            
              
                K
                D
              
              ⁢
              
                x
                .
              
            
          
        
      
    
    
      
        (
        32
        )
      
    
  




with Lag-stabilization

m{umlaut over (x)}+kx==uLS=−cx(t−τ).

Note that from PID the PD portions contribute to system stiffness and damping. What does the “I” or integral portion contribute to Negative damping. Since lag-stabilization was shown to phase shift x to {dot over (x)}, some prescribed amount of damping, then for the “integration” of x, the phase shift would be proportional to −{dot over (x)}. Recognizing this then





  
    α
    ⁢
    
      
        ∫
        o
        t
      
      ⁢
      
        x
        ⁢
        
          
        
        ⁢
        
          d
          τ
        
      
    
  
  =
  
    -
    
      
        x
        .
      
      .
    
  




Differentiating yields

αx=−{umlaut over (x)} whichimpliesthat {umlaut over (x)}+αx=0

for no net damping (neutral stability) and where





  α
  =
  
    
      ω
      2
    
    =
    
      
        
          K
          +
          
            K
            P
          
        
        M
      
      .
    
  




For the mass-spring-damper system with PID control (32) the stability boundary is at

−To{dot over (S)}i+{dot over (W)}=0

or

To{dot over (S)}i={dot over (W)}

where the averages can be removed for this linear system. Substituting the appropriate terms from (32) gives





  
    
      [
      
        
          (
          
            C
            +
            
              K
              D
            
          
          )
        
        ⁢
        
          x
          .
        
      
      ]
    
    ⁢
    
      x
      .
    
  
  =
  
    
      -
      
        [
        
          
            K
            I
          
          ⁢
          
            
              ∫
              o
              t
            
            ⁢
            
              x
              ⁢
              
                
              
              ⁢
              
                d
                τ
              
            
          
        
        ]
      
    
    ⁢
    
      x
      .
    
  




then rearranging as





  
    
      ∫
      o
      t
    
    ⁢
    
      x
      ⁢
      
        
      
      ⁢
      
        d
        τ
      
    
  
  =
  
    -
    
      
        (
        
          C
          +
          
            K
            D
          
        
        )
      
      
        K
        I
      
    
  




and differentiating both sides gives





  x
  =
  
    
      -
      
        
          (
          
            C
            +
            
              K
              D
            
          
          )
        
        
          K
          I
        
      
    
    ⁢
    
      x
       ̈
    
  




which results in

{umlaut over (x)}+ω2x=0 whichimpliesthat





  
    ω
    2
  
  =
  
    
      
        K
        I
      
      
        (
        
          C
          +
          
            K
            D
          
        
        )
      
    
    .
  




Next this result is further clarified with a conventional Routh-Hurwitz analysis.

First, convert the combined mass-spring-damper PID control system into the following equivalent third-order system




  
    
      M
      ⁢
      
        x
         ̈
      
    
    +
    
      
        (
        
          C
          +
          
            K
            D
          
        
        )
      
      ⁢
      
        x
        .
      
    
    +
    
      
        (
        
          K
          +
          
            K
            P
          
        
        )
      
      ⁢
      x
    
    +
    
      
        K
        I
      
      ⁢
      
        
          ∫
          o
          t
        
        ⁢
        
          x
          ⁢
          
            
          
          ⁢
          
            d
            τ
          
        
      
    
  
  =
  0




Next, invoke the following change of variable as





  y
  =
  
    
      ∫
      o
      t
    
    ⁢
    
      x
      ⁢
      
        
      
      ⁢
      
        d
        τ
      
    
  




gives





  
    
      y
      ...
    
    +
    
      
        
          (
          
            C
            +
            
              K
              D
            
          
          )
        
        M
      
      ⁢
      
        y
         ̈
      
    
    +
    
      
        
          (
          
            K
            +
            
              K
              P
            
          
          )
        
        M
      
      ⁢
      
        y
        .
      
    
    +
    
      
        
          K
          I
        
        M
      
      ⁢
      y
    
  
  =
  0




and transforming to the s-domain yields





  
    
      s
      3
    
    +
    
      
        
          (
          
            C
            +
            
              K
              D
            
          
          )
        
        M
      
      ⁢
      
        s
        2
      
    
    +
    
      
        
          (
          
            K
            +
            
              K
              P
            
          
          )
        
        M
      
      ⁢
      s
    
    +
    
      
        
          K
          I
        
        M
      
      ⁢
      y
    
  
  =
  0.




For the third-order system to be stable, the Routh-Hurwitz analysis results in the following necessary and sufficient conditions:

1. All the polynomial coefficients must be positive and
2.




  
    
      (
      
        
          (
          
            C
            +
            
              K
              D
            
          
          )
        
        M
      
      )
    
    ⁢
    
      (
      
        
          (
          
            K
            +
            
              K
              P
            
          
          )
        
        M
      
      )
    
  
  ≥
  
    
      
        K
        I
      
      M
    
    .
  




For the specific equality condition:





  
    
      (
      
        
          (
          
            C
            +
            
              K
              D
            
          
          )
        
        M
      
      )
    
    ⁢
    
      (
      
        
          (
          
            K
            +
            
              K
              P
            
          
          )
        
        M
      
      )
    
  
  =
  
    
      K
      I
    
    M
  




impliesthat





  
    
      K
      I
    
    
      (
      
        C
        +
        
          K
          D
        
      
      )
    
  
  =
  
    
      
        (
        
          K
          +
          
            K
            P
          
        
        )
      
      M
    
    =
    
      
        ω
        2
      
      .
    
  




This result is the marginal or neutral stability boundary which agrees with the previous lag stabilized and exergy/entropy control results:

α=ω2.

Also note that for the previous PID control numerical simulation results that case 1 is asymptotically stable, case 2 is neutrally stable, and case 3 is unstable which matches the analysis results.

EXAMPLE 3
This example is the design of a control law for a single degree of freedom nonlinear oscillator. The Duffing oscillator/Coulomb friction dynamic model (see FIG. 22) is defined as

M{umlaut over (x)}+C{dot over (x)}+CNL sign({dot over (x)})+Kx+KNLx3=u  (33)

where M, C, K, and u are the mass, damper, stiffness coefficients and external force input terms. The nonlinear stiffness and Coulomb friction coefficients are KNL and CNL, respectively. The PID controller is defined as





  
    
      
        u
        =
        
          
            
              -
              
                K
                p
              
            
            ⁢
            x
          
          -
          
            
              K
              I
            
            ⁢
            
              
                ∫
                0
                t
              
              ⁢
              
                x
                ⁢
                
                  
                
                ⁢
                
                  d
                  τ
                
              
            
          
          -
          
            
              K
              D
            
            ⁢
            
              x
              .
            
          
        
      
    
    
      
        (
        34
        )
      
    
  




where KP, KI, and KD are the proportional, integral and derivative controller gains.

Initially, the nonlinear Duffing oscillator is investigated as a neutrally stable, reversible conservative system or

M{umlaut over (x)}+Kx+KNLx3=−KPx 

subject to the initial condition x(0)=xo=1.0. Now apply exergy/entropy control design and the derivative of the Lyapunov function/Hamiltonian becomes





  
    
      
        
          V
          .
        
        =
        
          
            H
            .
          
          =
          
            
              
                W
                .
              
              -
              
                
                  T
                  o
                
                ⁢
                
                  
                    S
                    .
                  
                  i
                
              
            
            =
            
              
                
                  ∑
                  
                    j
                    =
                    1
                  
                  N
                
                ⁢
                
                  
                
                ⁢
                
                  
                    Q
                    j
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    j
                  
                
              
              -
              
                
                  ∑
                  
                    l
                    =
                    N
                  
                  
                    M
                    -
                    N
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    Q
                    l
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    l
                  
                
              
            
          
        
      
    
    
      
        (
        35
        )
      
    
  




which yields

To{dot over (S)}i=0

{dot over (W)}=0

(To{dot over (S)}rev)ave=(M{umlaut over (x)}·{dot over (x)}+(K+KP)x·{dot over (x)}+KNLx3·{dot over (x)})ave=0.

Numerical simulations are performed with the numerical values listed in Table 3. For this initial Case 1 the phase plane plot and the potential and kinetic energy rate plots are shown in FIGS. 23(a) and (b). This run demonstrates Corollary 1 and a stable orbit for the nonlinear system with offsetting potential and kinetic energy rates responses.







TABLE 3







Duffing oscillator/Coulomb friction 


model and PID control system gains















Case
KP
KI
KD
M
K
C
KNL
CNL


No.
(kg/s2)
(kg/s3)
(kg/s)
(kg)
(kg/s2)
(kg/s)
(N/m3)
(N)


















1
10.0
0.0
0.0
10.0
10.0
0.0
100.0
0.0


2
10.0
20.0
2.0
10.0
10.0
0.1
100.0
5.0


3
10.0
40.05
2.0
10.0
10.0
0.1
100.0
5.0


4
10.0
80.0
2.0
10.0
10.0
0.1
100.0
5.0









Next, consider the additional PID and Coulomb friction effects applied to the Duffing oscillator and partition into exergy generation and exergy dissipation terms. Now apply the exergy/entropy control design and the derivative of the Lyapunov function/Hamiltonian becomes




  
    
      
        
          V
          .
        
        =
        
          
            H
            .
          
          =
          
            
              
                W
                .
              
              -
              
                
                  T
                  o
                
                ⁢
                
                  
                    S
                    .
                  
                  i
                
              
            
            =
            
              
                
                  ∑
                  
                    j
                    =
                    1
                  
                  N
                
                ⁢
                
                  
                
                ⁢
                
                  
                    Q
                    j
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    j
                  
                
              
              -
              
                
                  ∑
                  
                    l
                    =
                    N
                  
                  
                    M
                    -
                    N
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    Q
                    l
                  
                  ⁢
                  
                    
                      q
                      .
                    
                    l
                  
                
              
            
          
        
      
    
    
      
        (
        36
        )
      
    
  




which yields





  
    
      
        
          
            T
            o
          
          ⁢
          
            
              S
              .
            
            i
          
        
        =
        
        ⁢
        
          
            
              (
              
                C
                +
                
                  K
                  D
                
              
              )
            
            ⁢
            
              
                x
                .
              
              ·
              
                x
                .
              
            
          
          +
          
            
              C
              NL
            
            ⁢
            
              
                sign
                ⁡
                
                  (
                  
                    x
                    .
                  
                  )
                
              
              ·
              
                x
                .
              
            
          
        
      
    
  
  
    
      
        
          W
          .
        
        =
        
        ⁢
        
          
            -
            
              K
              I
            
          
          ⁢
          
            
              ∫
              o
              t
            
            ⁢
            
              x
              ⁢
              
                
              
              ⁢
              
                
                  d
                  τ
                
                ·
                
                  x
                  .
                
              
            
          
        
      
    
  
  
    
      
        
          
            (
            
              
                T
                o
              
              ⁢
              
                
                  S
                  .
                
                rev
              
            
            )
          
          ave
        
        =
        
        ⁢
        
          
            
              (
              
                
                  M
                  ⁢
                  
                    
                      x
                       ̈
                    
                    ·
                    
                      x
                      .
                    
                  
                
                +
                
                  
                    (
                    
                      K
                      +
                      
                        K
                        p
                      
                    
                    )
                  
                  ⁢
                  
                    x
                    ·
                    
                      x
                      .
                    
                  
                
                +
                
                  
                    K
                    NL
                  
                  ⁢
                  
                    
                      x
                      3
                    
                    ·
                    
                      x
                      .
                    
                  
                
              
              )
            
            ave
          
          =
          0.
        
      
    
  




To determine the nonlinear stability boundary from the exergy/entropy control design

{dot over (V)}={dot over (H)}={dot over (W)}−To{dot over (S)}i 

which gives

({dot over (W)})ave=(To{dot over (S)}i)ave.

Substituting the actual terms yields the following:





  
    
      
        
          
            [
            
              
                -
                
                  K
                  I
                
              
              ⁢
              
                
                  ∫
                  o
                  t
                
                ⁢
                
                  x
                  ⁢
                  
                    
                  
                  ⁢
                  
                    
                      d
                      τ
                    
                    ·
                    
                      x
                      .
                    
                  
                
              
            
            ]
          
          ave
        
        =
        
          
            [
            
              
                
                  (
                  
                    C
                    +
                    
                      K
                      D
                    
                  
                  )
                
                ⁢
                
                  
                    x
                    .
                  
                  ·
                  
                    x
                    .
                  
                
              
              +
              
                
                  C
                  NL
                
                ⁢
                
                  
                    sign
                    (
                    
                      x
                      .
                    
                    )
                  
                  ·
                  
                    x
                    .
                  
                
              
            
            ]
          
          ave
        
      
    
    
      
        (
        37
        )
      
    
  




which is the nonlinear stability boundary. To best understand how the boundary is determined, concepts and analogies from electric AC power have been introduced earlier. Essentially, when the average powerin is equivalent to the average powerdissipated over a cycle, then the system is operating at the stability boundary. Later, in the exergy and exergy rate responses for the nonlinear system, one may observe that the area under the curves for the exergy rate generation and the exergy rate dissipation are equivalent and for the corresponding exergy responses the slopes will be equal and opposite. This helps to explain why PID control works well for nonlinear systems.

Numerical simulations are performed to demonstrate where the nonlinear stability boundary lies for the Duffing oscillator/Coulomb friction dynamic model subject to PID control. Three separate cases are conducted with the numerical values listed in Table 3. The nonlinear system is subject to an initial condition of x0=1.0. For case 2 the integral of position, position, velocity, and acceleration responses along with the exergy and exergy rate responses are plotted in FIGS. 24(a) and (b). For this case, the dissipative term is greater than the generative term. This is observed from the decaying system responses. In case 3 the system responses along with the exergy and exergy rate responses are shown in FIGS. 25(a) and (b). In this case, the average exergy slopes and integrated power areas for the dissipative and generative terms are equivalent which demonstrates (37). This results in system responses that do not decay, displaying constant nonlinear oscillatory behavior. In final case 4, the system responses along with the exergy and exergy rate responses are shown in FIGS. 26(a) and (b). In this case, the dissipative term is less than the generative term which results in a system response with increasing nonlinear oscillatory behavior. In conclusion, FIGS. 27(a) and (b) show the responses for the total exergy with respect to each case along with the phase plane plot for the nonlinear system. For case 3 the nonlinear stability boundary (or neutral stability) is characteristic of an average zero output for the total exergy response or validation of (37). For the phase plane plot, case 2 demonstrates an asymptotically stable decaying response, case 3 a neutrally stable orbital response, and case 4 an asymptotically unstable increasing orbit response.
All three cases for the PID control regulator Duffing oscillator/Coulomb friction dynamic system demonstrates the three subcases for Corollary 4: Given apriori (To{dot over (S)}i)ave>0 and ({dot over (W)})ave>0 then the nonlinear system, Duffing oscillator/Coulomb friction with PID control showed the following:

    
    
        4.1 Case 2 yielded (To{dot over (S)}i)ave>({dot over (W)})ave; asymptotic stability; damped stable nonlinear response and demonstration of Corollary 4.1.
        4.2 Case 3 yielded (To{dot over (S)}i)ave=({dot over (W)})ave; neutral stability; and demonstration of Corollary 4.2. This case is the nonlinear stability boundary where dissipation and generation terms cancel each other out on the average.
        4.3 Case 3 yielded (To{dot over (S)}i)ave<({dot over (W)})ave; increasingly unstable towards another orbit; and demonstration of Corollary 4.3.

In addition, the boundary condition can be used to identify different operating regions for systems that may need to be gain scheduled. For example, the integral gain, in general may be a function of many different parameters

KI=KI(xo,{dot over (x)}o,C,CNL,KD).

For the purpose of illustration the integral gain was investigated for different initial conditions, xo, while holding all the other possible parameters constant. The results are shown in FIG. 28. For several varying initial conditions for xo, at the stability boundary condition each KI was determined. The corresponding linear system resulted in constant KI values which is characteristic of a linear system. The gain scheduling is due to the nonlinear spring. It changes the limit cycle behavior as a function of store exergy. In the next example, exergy/entropy control design is applied to the tracking control problem.

    
    


EXAMPLE 4
The tracking controller design according to the invention begins with picking a Lyapunov function that is the error energy and the second variation of the Hamiltonian.




  
    
      
        V
        =
        
          
            1
            2
          
          ⁢
          
            δ
            2
          
          ⁢
          H
        
      
    
    
      
        (
        38
        )
      
    
  




which is positive definite. The time derivative is





  
    
      
        
          
            
              
                V
                .
              
              =
              
              ⁢
              
                
                  
                    d
                    
                      d
                      t
                    
                  
                  ⁢
                  
                    (
                    
                      
                        1
                        2
                      
                      ⁢
                      
                        δ
                        2
                      
                      ⁢
                      H
                    
                    )
                  
                
                =
                
                  δ
                  ⁢
                  
                    H
                    .
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                
                  
                    Δ
                    ⁢
                    
                      W
                      .
                    
                  
                  -
                  
                    
                      T
                      o
                    
                    ⁢
                    Δ
                    ⁢
                    
                      
                        S
                        .
                      
                      i
                    
                  
                
                =
                
                  
                    
                      ∑
                      
                        j
                        =
                        1
                      
                      N
                    
                    ⁢
                    
                      
                    
                    ⁢
                    
                      δ
                      ⁢
                      
                        
                      
                      ⁢
                      
                        Q
                        j
                      
                      ⁢
                      δ
                      ⁢
                      
                        
                          q
                          .
                        
                        j
                      
                    
                  
                  -
                  
                    
                      ∑
                      
                        l
                        =
                        N
                      
                      
                        M
                        -
                        N
                      
                    
                    ⁢
                    
                      
                    
                    ⁢
                    
                      δ
                      ⁢
                      
                        
                      
                      ⁢
                      
                        Q
                        l
                      
                      ⁢
                      δ
                      ⁢
                      
                        
                          q
                          .
                        
                        l
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        39
        )
      
    
  




that is subject to the following general necessary and sufficient conditions:

ToΔ{dot over (S)}i≧0 Positivesemi-definite, alwaystrue

Δ{dot over (W)}≧0 Positivesemi-definite; exergypumpedintothesystem.

Similar arguments can be made for the corollaries as compared to the earlier regulator problem cases.

The formation is identical to the regulator problem with the exception that for the nonlinear case additional cross-terms will need to be accommodated.
For the tracking controller, the variations are about a reference input trajectory where the Hamiltonian is viewed as

H=Hr+δH. 

This is best demonstrated by the following example. Initially, the mass-spring-damper (see FIG. 4) dynamic model (25) is revisited and a tracking PID control law is applied. For tracking control, first define the perturbation variables to be

x=xr+δx 

{dot over (x)}={dot over (x)}r+δ{dot over (x)}

{umlaut over (x)}={umlaut over (x)}r+δ{umlaut over (x)}

u=ur+δu 

with the reference control given as

ur=M{umlaut over (x)}r+C{dot over (x)}r+Kxr.

Then the PID tracking controller is





  
    δ
    ⁢
    
      
    
    ⁢
    u
  
  =
  
    
      -
      
        
          K
          P
        
        ⁡
        
          (
          
            δ
            ⁢
            
              
            
            ⁢
            x
          
          )
        
      
    
    -
    
      
        K
        I
      
      ⁢
      
        
          ∫
          o
          t
        
        ⁢
        
          
            (
            
              δ
              ⁢
              
                
              
              ⁢
              x
            
            )
          
          ⁢
          
            
          
          ⁢
          
            d
            τ
          
        
      
    
    -
    
      
        K
        D
      
      ⁡
      
        (
        
          δ
          ⁢
          
            x
            .
          
        
        )
      
    
  




where KP, KI, and KD are the proportional, integral, and derivative controller gains.

Now investigate the PID tracking control in terms of exergy dissipation and exergy generation. The derivative of the Lyapunov function/Hamiltonian becomes




  
    
      
        
          
            
              
                V
                .
              
              =
              
                
                  δ
                  ⁢
                  
                    H
                    .
                  
                
                =
                
                ⁢
                
                  
                    
                      Δ
                      ⁢
                      
                        W
                        .
                      
                    
                    -
                    
                      
                        T
                        o
                      
                      ⁢
                      Δ
                      ⁢
                      
                        
                          S
                          .
                        
                        i
                      
                    
                  
                  =
                  
                    
                      
                        ∑
                        
                          j
                          =
                          1
                        
                        N
                      
                      ⁢
                      
                        
                      
                      ⁢
                      
                        δ
                        ⁢
                        
                          
                        
                        ⁢
                        
                          Q
                          j
                        
                        ⁢
                        δ
                        ⁢
                        
                          
                            q
                            .
                          
                          j
                        
                      
                    
                    -
                    
                      
                        ∑
                        
                          l
                          =
                          N
                        
                        
                          M
                          -
                          N
                        
                      
                      ⁢
                      
                        
                      
                      ⁢
                      
                        δ
                        ⁢
                        
                          
                        
                        ⁢
                        
                          Q
                          l
                        
                        ⁢
                        δ
                        ⁢
                        
                          
                            q
                            .
                          
                          l
                        
                      
                    
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                
                  M
                  ⁢
                  
                    
                  
                  ⁢
                  δ
                  ⁢
                  
                    
                      x
                       ̈
                    
                    ·
                    δ
                  
                  ⁢
                  
                    x
                    .
                  
                
                +
                
                  
                    (
                    
                      K
                      +
                      
                        K
                        P
                      
                    
                    )
                  
                  ⁢
                  δ
                  ⁢
                  
                    
                  
                  ⁢
                  
                    x
                    ·
                    δ
                  
                  ⁢
                  
                    x
                    .
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                
                  
                    -
                    
                      (
                      
                        C
                        +
                        
                          K
                          D
                        
                      
                      )
                    
                  
                  ⁢
                  δ
                  ⁢
                  
                    
                      x
                      .
                    
                    ·
                    δ
                  
                  ⁢
                  
                    x
                    .
                  
                
                -
                
                  
                    K
                    I
                  
                  ⁢
                  
                    
                      ∫
                      o
                      t
                    
                    ⁢
                    
                      δ
                      ⁢
                      
                        
                      
                      ⁢
                      x
                      ⁢
                      
                        
                      
                      ⁢
                      
                        
                          d
                          τ
                        
                        ·
                        δ
                      
                      ⁢
                      
                        x
                        .
                      
                    
                  
                
              
            
          
        
      
    
    
      
        (
        40
        )
      
    
  




which yields





  
    
      
        
          
            
              
                
                  T
                  o
                
                ⁢
                Δ
                ⁢
                
                  
                
                ⁢
                
                  S
                  i
                
              
              =
              
              ⁢
              
                
                  (
                  
                    C
                    +
                    
                      K
                      D
                    
                  
                  )
                
                ⁢
                δ
                ⁢
                
                  
                    x
                    .
                  
                  ·
                  δ
                
                ⁢
                
                  x
                  .
                
              
            
          
        
        
          
            
              
                Δ
                ⁢
                
                  W
                  .
                
              
              =
              
              ⁢
              
                
                  -
                  
                    K
                    I
                  
                
                ⁢
                
                  
                    ∫
                    0
                    t
                  
                  ⁢
                  
                    δ
                    ⁢
                    
                      
                    
                    ⁢
                    x
                    ⁢
                    
                      
                    
                    ⁢
                    
                      
                        d
                        τ
                      
                      ·
                      δ
                    
                    ⁢
                    
                      x
                      .
                    
                  
                
              
            
          
        
        
          
            
              
                
                  [
                  
                    
                      T
                      o
                    
                    ⁢
                    Δ
                    ⁢
                    
                      
                        S
                        .
                      
                      rev
                    
                  
                  ]
                
                ave
              
              =
              
              ⁢
              
                
                  
                    [
                    
                      
                        M
                        ⁢
                        
                          
                        
                        ⁢
                        δ
                        ⁢
                        
                          
                            x
                             ̈
                          
                          ·
                          δ
                        
                        ⁢
                        
                          x
                          .
                        
                      
                      +
                      
                        
                          (
                          
                            K
                            +
                            
                              K
                              P
                            
                          
                          )
                        
                        ⁢
                        δ
                        ⁢
                        
                          
                        
                        ⁢
                        
                          x
                          ·
                          δ
                        
                        ⁢
                        
                          x
                          .
                        
                      
                    
                    ]
                  
                  ave
                
                =
                0.
              
            
          
        
      
    
    
      
        (
        41
        )
      
    
  




The first equation in (41), the derivative tracking control term along with the damping term are identified as dissipative terms. The second equation in (41), the integral tracking control term is identified as a generative term. As in the previous regulator case the same is true for the tracking case, in the final equation in (41), the proportional tracking control term is added to the system stiffness term and is identified as a reversible term.

Numerical simulations are performed for three separate PID tracking control cases with the numerical values listed in Table 4. The reference input signal is defined as

xr=Ain sin ωint. 

The mass-spring-damper system initial conditions are zero. For case 1 the tracking force input, δu, and tracking position and velocity errors are plotted in FIGS. 29(a) and (b). The variations for potential/kinetic energy rate, exergy, and exergy rate responses are shown in FIGS. 30(a) and (b). For case 1 passive stable tracking is observed from the decaying response and shows a predominately larger perturbation dissipative term with respect to the perturbation generative term. The initial transient tracking position and velocity errors eventually converge to the reference sinusoidal inputs. Case 1 corresponds with the tracking control equivalent to Corollary 4.1. For case 2, similar responses are given in FIGS. 31(a) and (b) for the tracking force input and tracking position and velocity errors. The variations for potential/kinetic energy rate, exergy, and exergy rate responses are shown in FIGS. 32(a) and (b). This case demonstrates the neutrally stable tracking boundary or where the dissipative terms are equivalent and cancel the generative terms which would correspond to the tracking control equivalent to Corollary 4.2. For case 3, similar responses are given in FIGS. 33(a) and (b) and 34(a) and (b). Case 3 shows responses that are growing exponentially, since the generative term is greater than the dissipative term. This final case would correspond to the tracking control equivalent to Corollary 4.3, respectively.








TABLE 4







Mass-spring-damper model and PID 


tracking control system gains















Case
KP
KI
KD
M/Mr
K/Kr
C/Cr
Ain
win


No.
(kg/s2) 
(kg/s3)
kg/s
(kg)
kg/s2)
(kg/s)
(m)
(rad/s)


















1
100.0
2.0
2.0
10.0/10.0 
10.0/10.0
0.1/0.1
1.0
1.0


2
100.0
23.1
2.0
10.0/10.0 
10.0/10.0
0.1/0.1
1.0
1.0


3
100.0
40.0
2.0
10.0/10.0 
10.0/10.0
0.1/0.1
1.0
1.0









The preceding examples can be repeated with similar success by substituting the generically or specifically described reactants and/or operating conditions of this invention for those used in the preceding examples.
Although the invention has been described in detail with particular reference to these preferred embodiments, other embodiments can achieve the same results. Variations and modifications of the present invention will be obvious to those skilled in the art and it is intended to cover in the appended claims all such modifications and equivalents. The entire disclosures of all references, applications, patents, and publications cited above are hereby incorporated by reference.

## Claims

1. A control system design method comprising the steps of:
representing a physical apparatus to be controlled as a Hamiltonian system, wherein the Hamiltonian system comprises a mass-spring-damper dynamic system;
determining elements of the Hamiltonian system representation which are power generators, power dissipators, and power storage devices;
analyzing stability and performance of the Hamiltonian system based on the results of the determining step and determining necessary and sufficient conditions for stability of the Hamiltonian system;
creating a stable control system based on the results of the analyzing step; and
employing the resulting control system to control the physical apparatus.

2. The method of claim 1 wherein the control system is nonlinear.

3. The method of claim 1 wherein in the analyzing step stability is determined in terms of power flow.

4. The method of claim 3 wherein in the analyzing step stability is determined in terms of whether the Hamiltonian system is moving toward or away from its minimum energy and maximum entropy state.

5. The method of claim 3 wherein in the analyzing step stability is determined via a Lyapunov derivative function.

6. The method of claim 5 where in the analyzing step stability is determined via a Lyapunov derivative function that is a decomposition and sum of exergy dissipation rate and exergy generation rate.

7. The method of claim 1 wherein the Hamiltonian system comprises a forced harmonic oscillator.

8. The method of claim 1 wherein the physical apparatus is selected from the group consisting of manufacturing devices, electronic devices, communications devices, transportation apparatuses, computer devices, computer networks, military apparatus, electrical power grids, Supervisory Control And Data Acquisition apparatus, and pipeline apparatus.

9. A control system for a physical apparatus, said control system comprising a stable control system based on results of analyzing stability and performance of the physical apparatus represented as a Hamiltonian system and based on the results of determining elements of the Hamiltonian system representation which are power generators, power dissipators, and power storage devices, and on the results of determining necessary and sufficient conditions for stability of the Hamiltonian system, wherein the Hamiltonian system comprises a mass-spring-damper dynamic system or a Duffing oscillator with Coulomb friction nonlinear damper system.

10. The control system of claim 9 wherein said control system is nonlinear.

11. The control system of claim 9 wherein stability is determined in terms of power flow.

12. The control system of claim 11 wherein stability is determined in terms of whether the Hamiltonian system is moving toward or away from its minimum energy and maximum entropy state.

13. The control system of claim 11 wherein stability is determined via a Lyapunov derivative function.

14. The control system of claim 13 where stability is determined via a Lyapunov derivative function that is a decomposition and sum of exergy dissipation rate and exergy generation rate.

15. The control system of claim 9 wherein the Hamiltonian system comprises a forced harmonic oscillator.

16. The control system of claim 9 wherein the physical apparatus is selected from the group consisting of manufacturing devices, electronic devices, communications devices, transportation apparatuses, computer devices, computer networks, military apparatus, electrical power grids, Supervisory Control And Data Acquisition apparatus, and pipeline apparatus.

17. A control system design method comprising the steps of:
representing a physical apparatus to be controlled as a Hamiltonian system, wherein the Hamiltonian system comprises a Duffing oscillator with Coulomb friction nonlinear damper system;
determining elements of the Hamiltonian system representation which are power generators, power dissipators, and power storage devices;
analyzing stability and performance of the Hamiltonian system based on the results of the determining step and determining necessary and sufficient conditions for stability of the Hamiltonian system;
creating a stable control system based on the results of the analyzing step; and
employing the resulting control system to control the physical apparatus.