System for intelligent control of a vehicle suspension based on soft computing

Information

  • Patent Grant
  • 6463371
  • Patent Number
    6,463,371
  • Date Filed
    Thursday, October 22, 1998
    25 years ago
  • Date Issued
    Tuesday, October 8, 2002
    21 years ago
Abstract
A reduced control system suitable for control of an active suspension system as a controlled plant is described. The reduced control system is configured to use a reduced sensor set for controlling the suspension without significant loss of control quality (accuracy) as compared to an optimal control system with an optimum sensor set. The control system calculates the information content provided by the reduced sensor set as compared to the information content provided by the optimum set. The control system also calculates the difference between the entropy production rate of the plant and the entropy production rate of the controller. A genetic optimizer is used to tune a fuzzy neural network in the reduced controller. A fitness function for the genetic optimizer provides optimum control accuracy in the reduced control system by minimizing the difference in entropy production while maximizing the sensor information content.
Description




BACKGROUND OF THE INVENTION




1. Field of the Invention




The disclosed invention relates generally to vehicle suspension systems, and more particularly to electronic control systems for active suspension systems.




2. Description of the Related Art




Feedback control systems are widely used to maintain the output of a dynamic system at a desired value in spite of external disturbance forces that would move the output away from the desired value. For example, a household furnace controlled by a thermostat is an example of a feedback control system. The thermostat continuously measures the air temperature of the house, and when the temperature falls below a desired minimum temperature, the thermostat turns the furnace on. When the furnace has warmed the air above the desired minimum temperature, then the thermostat turns the furnace off. The thermostat-furnace system maintains the household temperature at a constant value in spite of external disturbances such as a drop in the outside air temperature. Similar types of feedback control are used in many applications.




A central component in a feedback control system is a controlled object, otherwise known as a process “plant,” whose output variable is to be controlled. In the above example, the plant is the house, the output variable is the air temperature of the house, and the disturbance is the flow of heat through the walls of the house. The plant is controlled by a control system. In the above example, the control system is the thermostat in combination with the furnace. The thermostat-furnace system uses simple on-off feedback control to maintain the temperature of the house. In many control environments, such as motor shaft position or motor speed control systems, simple on-off feedback control is insufficient. More advanced control systems rely on combinations of proportional feedback control, integral feedback control, and derivative feedback control.




The PID control system is a linear control system that is based on a dynamic model of the plant. In classical control systems, a linear dynamic model is obtained in the form of dynamic equations, usually ordinary differential equations. The plant is assumed to be relatively linear, time invariant, and stable. However, many real-world plants are time varying, highly nonlinear, and unstable. For example, the dynamic model may contain parameters (e.g., masses, inductances, aerodynarmic coefficients, etc.) which are either poorly known or depend on a changing environment. If the parameter variation is small and the dynamic model is stable, then the PID controller may be sufficient. However, if the parameter variation is large, or if the dynamic model is unstable, then it is common to add adaptation or intelligent (AI) control to the PID control system.




AI control systems use an optimizer, typically a nonlinear optimizer, to program the operation of the PID controller and thereby improve the overall operation of the control system. The optimizers used in many AI control systems rely on a genetic algorithm. Using a set of inputs, and a fitness function, the genetic algorithm works in a manner similar to process of evolution to arrive at a solution which is, hopefully, optimal. The genetic algorithm generates sets of chromosomes (corresponding to possible solutions) and then sorts the chromosomes by evaluating each solution using the fitness function. The fitness function determines where each solution ranks on a fitness scale. Chromosomes which are more fit, are those chromosomes which correspond to solutions that rate high on the fitness scale. Chromosomes which are less fit, are those chromosomes which correspond to solutions that rate low on the fitness scale. Chromosomes that are more fit are kept (survive) and chromosomes that are less fit are discarded (die). New chromosomes are created to replace the discarded chromosomes. The new chromosomes are created by crossing pieces of existing chromosomes and by introducing mutations.




The PID controller has a linear transfer function and thus is based upon a linearized equation of motion for the plant. Prior art genetic algorithms used to program PID controllers typically use simple fitness functions and thus do not solve the problem of poor controllability typically seen in linearization models. As is the case with most optimizers, the success or failure of the optimization often ultimately depends on the selection of the performance (fitness) function.




Evaluating the motion characteristics of a nonlinear plant is often difficult, in part due to the lack of a general analysis method. Conventionally, when controlling a plant with nonlinear motion characteristics, it is common to find certain equilibrium points of the plant and the motion characteristics of the plant are linearized in a vicinity near an equilibrium point. Control is then based on evaluating the pseudo (linearized) motion characteristics near the equilibrium point. This technique works poorly, if at all, for plants described by models that are unstable or dissipative.




SUMMARY OF THE INVENTION




The present invention solves these and other problems by providing a new AI control system that allows a reduced number of sensors to be used without a significant loss in control accuracy. The new AI control system is self-organizing and uses a fitness (performance) function that are based on the physical laws of minimum entropy and maximum sensor information. The self-organizing control system may be used to control complex plants described by nonlinear, unstable, dissipative models. The reduced control system is configured to use smart simulation techniques for controlling the plant despite the reduction in the number of sensor number without significant loss of control quality (accuracy) as compared to an optimal control system. In one embodiment, the reduced control system comprises a neural network that is trained by a genetic analyzer. The genetic analyzer uses a fitness function that maximizes information while minimizing entropy production.




In one embodiment, the reduced control system is applied to an internal combustion engine to provide control without the use of extra sensors, such as, for example, an oxygen sensor. The reduced control system develops a reduced control signal from a reduced sensor set. The reduced control system is trained by a genetic analyzer that uses a control signal developed by an optimized control system. The optimized control system provides an optimum control signal based on data obtained from temperature sensors, air-flow sensors, and an oxygen sensor. In an off-line learning mode, the optimum control signal is subtracted from a reduced control signal (developed by the reduced control system) and provided to an information calculator. The information calculator provides an information criteria to the genetic analyzer. Data from the reduced sensor set is also provided to an entropy model, which calculates a physical criteria based on entropy. The physical criteria is also provided to the genetic analyzer. The genetic analyzer uses both the information criteria and the physical criteria to develop a training signal for the reduced control system.




In one embodiment, a reduced control system is applied to a vehicle suspension to provide control of the suspension system using data from a reduced number of sensors. The reduced control system develops a reduced control signal from a reduced sensor set. The reduced control system is trained by a genetic analyzer that uses a control signal developed by an optimized control system. The optimized control system provides an optimum control signal based on data obtained from a plurality of angle and position sensors. In an off-line learning mode, the optimum control signal is subtracted from a reduced control signal (developed by the reduced control system) and provided to an information calculator. In one embodiment, the reduced control system uses a vertical accelerometer mounted near the center of the vehicle. The information calculator provides an information criteria to the genetic analyzer. Data from the reduced sensor set is also provided to an entropy model, which calculates a physical criteria based on entropy. The physical criteria is also provided to the genetic analyzer. The genetic analyzer uses both the information criteria and the physical criteria to develop a training signal for the reduced control system.




In one embodiment, the invention includes a method for controlling a nonlinear object (a plant) by obtaining an entropy production difference between a time differentiation (dS


u


/dt) of the entropy of the plant and a time differentiation (dS


c


/dt) of the entropy provided to the plant from a controller. A genetic algorithm that uses the entropy production difference as a fitness (performance) function evolves a control rule for a low-level controller, such as a PID controller. The nonlinear stability characteristics of the plant are evaluated using a Lyapunov function. The evolved control rule may be corrected using further evolutions using an information function that compares the information available from an optimum sensor system with the information available from a reduced sensor system. The genetic analyzer minimizes entropy and maximizes sensor information content.




In some embodiments, the control method may also include evolving a control rule relative to a variable of the controller by means of a genetic algorithm. The genetic algorithm uses a fitness function based on a difference between a time differentiation of the entropy of the plant (dS


u


/dt) and a time differentiation (dS


c


/dt) of the entropy provided to the plant. The variable may be corrected by using the evolved control rule.




In another embodiment, the invention comprises an AI control system adapted to control a nonlinear plant. The AI control system includes a simulator configured to use a thermodynamic model of a nonlinear equation of motion for the plant. The thermodynamic model is based on a Lyapunov function (V), and the simulator uses the function V to analyze control for a state stability of the plant. The AI control system calculates an entropy production difference between a time differentiation of the entropy of said plant (dS


u


/dt) and a time differentiation (dS


c


/dt) of the entropy provided to the plant by a low-level controller that controls the plant. The entropy production difference is used by a genetic algorithm to obtain an adaptation function in which the entropy production difference is minimized. The genetic algorithm provides a teaching signal to a fuzzy logic classifier that determines a fuzzy rule by using a learning process. The fuzzy logic controller is also configured to form a control rule that sets a control variable of the low-level controller.




In one embodiment, the low-level controller is a linear controller such as a PID controller. The learning processes may be implemented by a fuzzy neural network configured to form a look-up table for the fuzzy rule.




In yet another embodiment, the invention comprises a new physical measure of control quality based on minimum production entropy and using this measure for a fitness function of genetic algorithm in optimal control system design. This method provides a local entropy feedback loop in the control system. The entropy feedback loop provides for optimal control structure design by relating stability of the plant (using a Lyapunov function) and controllability of the plant (based on production entropy of the control system). The control system is applicable to all control systems, including, for example, control systems for mechanical systems, bio-mechanical systems, robotics, electromechanical systems, etc.











BRIEF DESCRIPTION OF THE DRAWINGS




The advantages and features of the disclosed invention will readily be appreciated by persons skilled in the art from the following detailed description when read in conjunction with the drawings listed below.





FIG. 1

is a block diagram showing an example of an AI control method in the prior art.





FIG. 2

is a block diagram showing an embodiment of an AI control method in accordance with one aspect of the present invention.





FIG. 3A

is a block diagram of an optimal control system





FIG. 3B

is a block diagram of a reduced control system.





FIG. 4A

is an overview block diagram of a reduced control system that uses entropy-based soft computing.





FIG. 4B

is a detailed block diagram of a reduced control system that uses entropy-based soft computing.





FIG. 5A

is a system block diagram of a fuzzy logic controller.





FIG. 5B

is a computing block diagram of the fuzzy logic controller shown in FIG.


5


A.





FIG. 6

is a diagram of an internal combustion piston engine with sensors.





FIG. 7

is a block diagram of a control system for controlling the internal combustion piston engine shown in FIG.


6


.





FIG. 8

is a schematic diagram of one half of an automobile suspension system.





FIG. 9

is a block diagram of a control system for controlling the automobile suspension system shown in FIG.


8


.











DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT




A feedback control system is commonly used to control an output variable of a process or plant in the face of some disturbance. Linear feedback control systems typically use combinations of proportional feedback control, integral feedback control, and derivative feedback control. Feedback that is the sum of proportional plus integral plus derivative feedback is often referred to as PID control. The Laplace transform of an output u(s) of a PBD controller is given by:










u


(
s
)


=



G


(
s
)




e


(
s
)



=


[


k
1

+


k
2

s

+


k
3


s


]



e


(
s
)








(
1
)













In the above equation, G(s) is the transfer function of the PID controller, e(s) is the controller input, u(s) is the controller output, k


1


is the coefficient for proportional feedback, k


2


is the coefficient for integral feedback, and k


3


is the coefficient for derivative feedback. The coefficients k


i


may be represented by a coefficient vector K, where K=[k


1


, k


2


, k


3


]. The vector K is commonly called a Coefficient Gain Schedule (CGS). The values of the coefficients K used in the linear PID control system are based on a dynamic model of the plant. When the plant is unstable, nonlinear, and/or time-variant, then the coefficients in K are often controlled by an AI control system.





FIG. 1

shows a typical prior art AI control system


100


. An input y(t) of the control system


100


is provided to a plus input of an adder


104


and an output x(t) of a plant


110


is provided to a minus input of the adder


104


. An output of the adder


104


is provided as an error signal e(t) to an error signal input of a PID controller


106


. An output u(t) of the PBD controller


106


is provided to a first input of an adder


108


. A disturbance m(t) is provided to a second input of the adder


108


. An output u*(t) of the adder


108


is provided to an input of the plant


110


. The plant


110


has a transfer function H(s) and an output x(t), where x(t)X(s) (where the symbol denotes the Laplace transform) and X(s)=H(s)u*(s). An output of the genetic algorithm


116


is provided to an input of a Fuzzy logic Neural Network (FNN)


118


and an output of the fuzzy neural network


118


is provided to a Fuzzy Controller (FC)


120


. An output of the fuzzy controller


120


is a set of coefficients K, which are provided to a coefficient input of the PID controller


106


.




The error signal e(t) provided to the PID controller


106


is the difference between the desired plant output value y(t) and the actual plant output value x(t). The PID controller


106


is designed to minimize the error represented by e(t) (the error being the difference between the desired and actual output signal signals). The PID controller


106


minimizes the error e(t) by generating an output signal u(t) which will move the output signal x(t) from the plant


110


closer to the desired value. The genetic algorithm


116


, fuzzy neural network


118


, and fuzzy controller


120


monitor the error signal e(t) and modify the gain schedule K of the PID controller


106


in order to improve the operation of the PID controller


106


.




The PID controller


106


constitutes a reverse model relative to the plant


110


. The genetic algorithm


116


evolves an output signal α based on a performance functions ƒ. Plural candidates for α are produced and these candidates are paired according to which plural chromosomes (parents) are produced. The chromosomes are evaluated and sorted from best to worst by using the performance function ƒ. After the evaluation for all parent chromosomes, good offspring chromosomes are selected from among the plural parent chromosomes, and some offspring chromosomes are randomly selected. The selected chromosomes are crossed so as to produce the parent chromosomes for the next generation. Mutation may also be provided. The second-generation parent chromosomes are also evaluated (sorted) and go through the same evolutionary process to produce the next-generation (i.e., third-generation) chromosomes. This evolutionary process is continued until it reaches a predetermined generation or the evaluation function ƒ finds a chromosome with a certain value. The outputs of the genetic algorithm are the chromosomes of the last generation. These chromosomes become input information α provided to the fuzzy neural network


118


.




In the fuzzy neural network


118


, a fuzzy rule to be used in the fuzzy controller


120


is selected from a set of rules. The selected rule is determined based on the input information α from the genetic algorithm


116


. Using the selected rule, the fuzzy controller


120


generates a gain schedule K for the PID controller


106


. The vector coefficient gain schedule K is provided to the PID controller


106


and thus adjusts the operation of the PID controller


106


so that the PID controller


106


is better able to minimize the error signal e(t).




Although the AI controller


100


is advantageous for accurate control in regions near linearized equilibrium points, the accuracy deteriorates in regions away from the linearized equilibrium points. Moreover, the AI controller


100


is typically slow or even unable to catch up with changes in the environment surrounding the plant


110


. The PID controller


106


has a linear transfer function G(s) and thus is based upon a linearized equation of motion for the plant


110


. Since the evaluation function ƒ used in the genetic algorithm


116


is only based on the information related to the input e(t) of the linear PID controller


106


, the controller


100


does not solve the problem of poor controllability typically seen in linearization models. Furthermore, the output results, both in the gain schedule K and the output x(t) often fluctuate greatly, depending on the nature of the performance function ƒ used in the genetic algorithm


116


. The genetic algorithm


116


is a nonlinear optimizer that optimizes the performance function ƒ. As is the case with most optimizers, the success or failure of the optimization often ultimately depends on the selection of the performance function ƒ.




The present invention solves these and other problems by providing a new AI control system


200


shown in FIG.


2


. Unlike prior AI control systems, the control system


200


is self-organizing and uses a new performance function ƒ, which is based on the physical law of minimum entropy. An input y(t) of the control system


200


is provided to a plus input of an adder


204


and an output x(t) of a plant


210


is provided to a minus input of the adder


204


. An output of the adder


204


is provided as an error signal e(t) to an error signal input of a PID controller


206


and to an input of a fuzzy controller


220


. An output u(t) of the PID controller


206


is provided to a first input of an adder


208


and to a first input of an entropy calculator (EC)


214


. A disturbance m(t) is provided to a second input of the adder


208


. An output u*(t) of the adder


208


is provided to an input of the plant


210


. The plant


210


has a transfer function H(s) and an output x(t), such that X(s)=H(s)u*(s), where x(t)X(s). The output x(t) is provided to a second input of the entropy calculator


214


and to the minus input of the adder


204


. An output of the entropy calculator


214


is provided to an input of a genetic algorithm


216


and an output of the genetic algorithm


216


is provided to an input of a Fuzzy logic Neural Network (FNN)


218


. An output of the fuzzy neural network


218


is provided to a rules selector input


222


of the fuzzy controller


220


. A Coefficient Gain Schedule (CGS) output


212


of the fuzzy controller


222


is provided to a gain schedule input of the PID


206


.




The combination of the genetic algorithm


216


and the entropy calculator


214


comprises a Simulation System of Control Quality


215


. The combination of the fuzzy neural network


218


and the fuzzy controller


220


comprises a Fuzzy Logic Classifier System FLCS


219


. The combination of the plant


210


and the adder


208


comprises a disturbed plant model


213


. The disturbed plant signal u*(t)=u(t)+m(t), and the disturbance m(t) are typically unobservable.




The error signal e(t) provided to the PID controller


206


is the difference between the desired plant output value y(t) and the actual plant output value x(t). The PID controller


206


is designed to minimize the error represented by e(t). The PID controller


206


minimizes the error e(t) by generating an output signal u(t) which will move the output signal x(t) from the plant


210


closer to the desired value. The fuzzy controller


220


monitors the error signal e(t) and modifies the gain schedule K of the PID controller


206


according to a fuzzy control rule selected by the fuzzy neural network


218


.




The genetic algorithm


216


provides a teaching signal K


T


to the fuzzy neural network


218


. The teaching signal K


T


is a global optimum solution of a coefficient gain schedule K generated by the genetic algorithm


216


.




The PID controller


206


constitutes a reverse model relative to the plant


210


. The genetic algorithm


216


evolves an output signal α based on a performance function ƒ. Plural candidates for α are produced and these candidates are paired by which plural chromosomes (parents) are produced. The chromosomes are evaluated and sorted from best to worst by using the performance function ƒ. After the evaluation for all parent chromosomes, good offspring chromosomes are selected from among the plural parent chromosomes, and some offspring chromosomes are randomly selected. The selected chromosomes are crossed so as to produce the parent chromosomes for the next generation. Mutation is also employed. The second-generation parent chromosomes are also evaluated (sorted) and go through the same evolutionary process to produce the next-generation (i.e., third-generation) chromosomes. This evolutionary process is continued until it reaches a predetermined generation or the evaluation function ƒ finds a chromosome with a certain value. Then, a component from a chromosome of the last generation becomes a last output, i.e., input information α provided to the fuzzy neural network


218


.




In the fuzzy neural network


218


, a fuzzy rule to be used in the fuzzy controller


220


is selected from a set of rules. The selected rule is determined based on the input information α from the genetic algorithm


216


. Using the selected rule, the fuzzy controller


220


generates a gain schedule K for the PID controller


206


. This is provided to the PID controller


206


and thus adjusts the operation of the PID controller


206


so that the PID controller


206


is better able to minimize the error signal e(t).




The fitness function ƒ for the genetic algorithm is given by









f
=

min








S



t







(
2
)













where












S



t


=

(





S
c




t


-




S
u




t



)





(
3
)













The quantity dS


u


/dt represents the rate of entropy production in the output x(t) of the plant


210


. The quantity dS


c


/dt represents the rate of entropy production in the output u(t) of the PID controller


206


.




Entropy is a concept that originated in physics to characterize the heat, or disorder, of a system. It can also be used to provide a measure of the uncertainty of a collection of events, or, for a random variable, a distribution of probabilities. The entropy function provides a measure of the lack of information in the probability distribution. To illustrate, assume that p(x) represents a probabilistic description of the known state of a parameter, that p(x) is the probability that the parameter is equal to z. If p(x) is uniform, then the parameter p is equally likely to hold any value, and an observer will know little about the parameter p. In this case, the entropy function is at its maximum. However, if one of the elements of p(z) occurs with a probability of one, then an observer will know the parameter p exactly and have complete information about p. In this case, the entropy of p(x) is at its minimum possible value. Thus, by providing a measure of uniformity, the entropy function allows quantification of the information on a probability distribution.




It is possible to apply these entropy concepts to parameter recovery by maximizing the entropy measure of a distribution of probabilities while constraining the probabilities so that they satisfy a statistical model given measured moments or data. Though this optimization, the distribution that has the least possible information that is consistent with the data may be found. In a sense, one is translating all of the information in the data into the form of a probability distribution. Thus, the resultant probability distribution contains only the information in the data without imposing additional structure. In general, entropy techniques are used to formulate the parameters to be recovered in terms of probability distributions and to describe the data as constraints for the optimization. Using entropy formulations, it is possible to perform a wide range of estimations, address ill-posed problems, and combine information from varied sources without having to impose strong distributional assumptions.




Entropy-based control of the plant in

FIG. 2

is based on obtaining the difference between a time differentiation (dS


u


/dt) of the entropy of the plant and a time differentiation (dS


c


/dt) of the entropy provided to the plant from a low-level controller that controls the plant


210


, and then evolving a control rule using a genetic algorithm. The time derivative of the entropy is called the entropy production rate. The genetic algorithm uses the difference between the entropy production rate of the plant (dS


u


/dt) and the entropy production rate of the low-level controller (dS


c


/dt) as a performance function. Nonlinear operation characteristics of the physical plant


210


are calculated by using a Lyapunov function




The dynamic stability properties of the plant


210


near an equilibrium point can be determined by use of Lyapunov functions. Let V(x) be a continuously differentiable scalar function defined in a domain D⊂R


n


that contains the origin. The function V(x) is said to be positive definite if V(0)=0 and V(x)>0 for x≠0. The function V(x) is said to be positive semidefinite if V(x)≧0 for all x. A function V(x) is said to be negative definite or negative semidefinite if −V(x) is positive definite or positive semidefinite, respectively. The derivative of V along the trajectories {dot over (x)}=ƒ(x) is given by:











V
.



(
x
)


=





i
=
1

n










V




x
i






x
.

i



=




V



x




f


(
x
)








(
4
)













where ∂V/∂x is a row vector whose ith component is ∂V/∂x


i


and the components of the n-dimensional vector ƒ(x) are locally Lipschitz functions of x, defined for all x in the domain D. The Lyapunov stability theorem states that the origin is stable if there is a continuously differentiable positive definite function V(x) so that {dot over (V)}(x) is negative definite. A function V(x) satisfying the conditions for stability is called a Lyapunov function.




Calculation of the Lyapunov dynamic stability and entropy production for a closed nonlinear mechanical system is demonstrated by using the Holmes-Rand (Duffing-Van der Pol) nonlinear oscillator as an example. The Holmes-Rand oscillator is described by the equation:






{umlaut over (x)}+(α+βx


2


){dot over (x)}−γx+x


3


=0  (5)






where α, β, and γ are constant parameters. A Lyapunov function for the Holmes-Rand oscillator is given by:










V
=



1
2




x
.

2


+

U


(
x
)




,






where





U

=



1
4



x
4


-


1
2


γ






x
2








(
6
)













Entropy production d


i


S/dt for the Holmes-Rand oscillator is given by the equation:













i


S



t


=


(

α
+

β






x
2



)




x
.

2






(
7
)













Equation 5 can be rewritten as:











x
¨

+


(

α
+

β






x
2



)



x
.


+



U



x



=
0




(
8
)













After multiplying both sides of the above equation by {dot over (x)}, then dV/dt can be calculated as:












V



t


=




x
¨







x
.


+




U



x




x
.



=


-





1
T







i


S



t








(
9
)













where T is a normalizing factor.




An interrelation between a Lyapunov function and the entropy production in an open dynamic system can be established by assuming a Lyapunov function of the form









V
=


1
2






i
=
1

n







(


q
i
2

+

S
2


)







(
10
)













where S=S


u


−S


c


, and the q


i


are generalized coordinates. It is possible to introduce the entropy function S in the Lyapunov function V because entropy S is also a scalar function of time. Differentiation of V with respect to time gives:












V



t


=





i
=
1

6








q
i




q
.

i



+

S


S
.







(
11
)













In this case, {dot over (q)}


i





i


(q


i


, τ, t), S=S


u


−S


c


, {dot over (S)}={dot over (S)}


u


−{dot over (S)}


c


and thus:












V



t


=





i
=
1

6








q
i




φ
i



(


q
i

,
τ
,
t

)




+


(


S
u

-

S
c


)



(





S
u




t


-




S
c




t



)







(
12
)













A special case occurs when β=0 and the Holmes-Rand oscillator reduces to a force-free Duffing, oscillator, wherein:













i


S



t


=


-
α








x
.

2







(

Duffing





oscillator

)






(
13
)













A Van der Pol oscillator is described by the equation:






{umlaut over (x)}+(x


2


−1){dot over (x)}+x=0  (14)






and the entropy production is given by:













i


S



t


=


1
T



(


x
2

-
1

)




x
.

2







(

Van





der





Pol





oscillator

)






(
15
)













For a micro-mobile robot in fluid, a mechanical model is given by:













m
1




x
¨

1


+


C
d



ρ
2



A
1



|


x
.

1

|



x
.

1

+


K
1



(


x
1

-

x
0

-


l
1



θ
0



)


-


K
2



(


x
2

-

x
1

-


l
2



θ
1



)




=
0




(
16
)























m
1




x
¨

2


+


C
d



ρ
2



A
2



|


x
.

2

|



x
.

2

+


K
2



(


x
2

-

x
1

-


l
2



θ
1



)


-


K
3



(


x
3

-

x
2

-


l
3



θ
2



)




=
0




(
17
)











m
3




x
¨

3


+


C
d



ρ
2



A
3



|


x
.

3

|



x
.

3

+


K
3



(


x
3

-

x
2

-


l
3



θ
2



)




=
0









(
18
)













where:










θ

n
+
1


=



-

1
2




θ
n


+


3
2







1

l

n
+
1





(


x

n
+
1


-

x
n


)







(
19
)













Values for a particular micro-mobile robot are given in Table 1 below.














TABLE 1









Item




Value




Units











m


1






1.6 × 10


−7






kg






m


2






1.4 × 10


−6






kg






m


3






2.4 × 10


−6






kg






l


1






2.0 × 10


−3






m






l


1






4.0 × 10


−3






m






l


3






4.0 × 10


−3






m






K


1






61.1




N/m






K


2






13.7




N/m






K


3






23.5




N/m






A


1






4.0 × 10


−6






m


2








A


2






2.4 × 10


−6






m


2








A


3






4.0 × 10


−6






m


2








C


d






 1.12











ρ




1000   




Kg/m


3
















Entropy production for the micro-mobile robot is given by the equation:













S
i




t


=




n
=
1

3








C
d



ρ
2



A
n



&LeftBracketingBar;


x
.

n

&RightBracketingBar;




x
.

n
2







(
20
)













and the Lyapunov Function is given by:









V
=





i
=
1

3








m
i




x
.

i
2



ρ
2



+




i
=
1

3










K
i



(


x
i

-

x

i
-
1


-


l
i



θ

i
-
1




)


2

2


+


S
2

2






(
21
)













where S=S


i


−S


c


and S


c


is the entropy of a controller with torque τ.




The necessary and sufficient conditions for Lyapunov stability of a plant is given by the relationship:













i




q
i




φ
i



(


q
i

,
τ
,
t

)




<


(


S
i

-


S
.

c


)



(





S
c




t


-




S
i




t



)



,





S
c




t


>




S
i




t







(
22
)













According to the above equation, stability of a plant can be achieved with “negentropy” S


c


(in the terminology used by Brillouin) where a change of negentropy dS


c


/dt in the control system


206


is subtracted from a change of entropy dS


i


/dt in the motion of the plant


210


.




The robust AI control system


200


provides improved control of mechanical systems in stochastic environments (e.g., active vibration control), intelligent robotics and electro-mechanical systems (e.g., mobile robot navigation, manipulators, collective mobile robot control), bio-mechanical systems (e.g., power assist systems, control of artificial replaced organs in medical systems as artificial lung ventilation), micro electromechanical systems (e.g., micro robots in fluids), etc.




The genetic algorithm realizes the search of optimal controllers with a simple structure using the principle of minimum entropy production. The fuzzy neural network controller offers a more flexible structure of controllers with a smaller torque, and the learning process produces less entropy. The fuzzy neural network controller gives a more flexible structure to controllers with smaller torque and the learning process produces less entropy than a genetic analyzer alone. Thus, an instinct mechanism produces less entropy than an intuition mechanism. However, necessary time for achieving an optimal control with learning process on fuzzy neural network (instinct) is larger than with the global search on genetic algorithm (intuition).




Realization of coordinated action between the look-up tables of the fuzzy controller


220


is accomplished by the genetic algorithm and the fuzzy neural network. In particular, the structure


200


provides a multimode fuzzy controller coupled with a linear or nonlinear neural network


218


. The control system


200


is a realization of a self-organizing AI control system with intuition and instinct. In the adaptive controller


200


, the feedback gains of the PID controller


210


are changed according to the quantum fuzzy logic, and approximate reasoning is provided by the use of nonlinear dynamic motion equations.




The fuzzy tuning rules for the gains k


i


are shaped by the learning system in the fuzzy neural network


218


with acceleration of fuzzy rules on the basis of global inputs provided by the genetic algorithm


216


. The control system


200


is thus a hierarchical, two-level control system that is intelligent “in small.” The lower (execution) level is provided by a traditional PID controller


206


, and the upper (coordination) level is provided by a KB (with fuzzy inference module in the form of production rules with different model of fuzzy implication) and fuzzification and de-fuzzification components, respectively.




The genetic algorithm


216


simulates an intuition mechanism of choosing the optimal structure of the PID controller


206


by using the fitness function, which is the measure of the entropy production, and the evolution function, which in this case is entropy.




Reduced Sensor Control Systems




In one embodiment, the entropy-based control system is applied to a reduced control system wherein the number of sensors has been reduced from that of an optimal system.





FIG. 3A

shows an optimal control system x


302


controlling a plant


304


. The optimal control system x produces an output signal x. A group of m sensors collect state information from the plant


304


and provide the state information as feedback to the optimal control system x


302


. The control system x


302


is optimal not in the sense that it is perfect or best, but rather, the control system x


302


is optimal in the sense that it provides a desired control accuracy for controlling the output of the plant


304


.





FIG. 3B

shows a reduced control system y


312


controlling the plant


304


. The reduced control system y


312


provides an output signal y. A group of n sensors (where n<m) collect state information from the plant


304


and provide the state information as feedback to the reduced control system x


312


. The reduced control system


312


is reduced because it receives information from fewer sensors than the optimal control system


302


(n is less than m).




The reduced control system


312


is configured to use smart simulation techniques for controlling the plant


304


despite the reduction in the number of sensor number without significant loss of control quality (accuracy) as compared to the optimal control system


302


. The plant


304


may be either stable or unstable, and the control system


312


can use either a complete model of the plant


304


or a partial model of the plant


304


.




The reduced control system


312


provides reliability, sufficient accuracy, robustness of control, stability of itself and the plant


304


, and lower cost (due, in part, to the reduced number of sensors needed for control). The reduced control system


312


is useful for robust intelligent control systems for plants such as, for example, suspension system, engines, helicopters, ships, robotics, micro electro-mechanical systems, etc.




In some embodiments, reduction in the number of sensors may be used, for example, to integrate sensors as intelligent systems (e.g., “sensor-actuator-microprocessor” with the preliminary data processing), and for simulation of robust intelligent control system with a reduced number of sensors.




The amount of information in the reduced output signal y (with n sensors) is preferably close to the amount of information in the optimal output signal. In other words, the reduced (variable) output y is similar to the optimal output x, such that (x≈y). The estimated accuracy ε of the output y as related to x is M[(x−y )


2


]≦ε


2


.




In some embodiments, the reduced control system


312


is adapted to minimize the entropy production in the control system


312


and in the plant


304


.




The structure of a reduced control system is shown in

FIGS. 4A and 4B

.

FIG. 4A

is a block diagram showing a reduced control system


480


and an optimal control system


420


. The optimal control system


420


, together with an optimizer


440


and a sensor compensator


460


are used to teach the reduced control system


480


. In

FIG. 4A

, a desired signal (representing the desired output) is provided to an input of the optimal control system


420


and to an input of the reduced control system


480


. The optimal control system


420


, having m sensors, provides an output sensor signal x


b


and an optimal control signal x


a


. A reduced control system


480


provides an output sensor signal y


b


and a reduced control signal y


a


. The signals x


b


and y


b


includes data from k sensors where k≦m−n. The k sensors typically being the sensors that are not common between the sensor systems


422


and


482


. The signals x


b


and y


b


are provided to first and second inputs of a subtractor


491


. An output of the subtractor


491


is a signal ε


b


, where ε


b


=x


b


−y


b


. The signal ε


b


is provided to a sensor input of a sensor compensator


460


. The signals x


a


and y


a


are provided to first and second inputs of a subtractor


490


. An output of the subtractor


490


is a signal ε


a


where ε


a


=x


a


−y


a


. The signal ε


a


is provided to a control signal input of the sensor compensator


460


. A control information output of the sensor compensator


460


is provided to a control information input of the optimizer


440


. A sensor information output of the sensor compensator


460


is provided to a sensor information input of the optimizer


440


. A sensor signal


483


from the reduced control system


480


is also provided to an input of the optimizer


440


. An output of the optimizer


440


provides a teaching signal


443


to an input of the reduced control system


480


.




In the description that follows, off-line mode typically refers to a calibration mode, wherein the control object


428


(and the control object


488


) is run with an optimal set of set of m sensors. In one embodiment, the off-line mode is run at the factory or repair center where the additional sensors (i.e., the sensors that are in the m set but not in the n set) are used to train the FNN


1




426


and the FNN


2




486


. The online mode typically refers to an operational mode (e.g., normal mode) where the system is run with only the n set of sensors.





FIG. 4B

is a block diagram showing details of the blocks in FIG.


4


A. In

FIG. 4B

, the output signal x


k


is provided by an output of a sensor set m


422


having m sensors given by m=k+k


1


. The information from the sensor system m


422


is a signal (or group of signals) having optimal information content I


1


. In other words, the information I


1


is the information from the complete set of m sensors in the sensor system


422


. The output signal x


a


is provided by an output of a control unit


425


. The control signal x


a


is provided to an input of a control object


428


. An output of the control object


428


is provided to an input of the sensor system


422


. Information I


k


from the set k of sensors is provided to an online learning input of a fuzzy neural network (FNN


1


)


426


and to an input of a first Genetic Algorithm (GA


1


)


427


. Information I


k1


from the set of sensors k


1


in the sensor system


422


is provided to an input of a control object model


424


. An off-line tuning signal output from the algorithm GA


1




427


is provided to an off-line tuning signal input of the FNN


1




426


. A control output from the FNN


426


is the control signal x


a


, which is provided to a control input of the control object


428


. The control object model


424


and the FNN


426


together comprise an optimal fuzzy control unit


425


.




Also in

FIG. 4B

, the sensor compensator


460


includes a multiplier


462


, a multiplier


466


, an information calculator


468


, and an information calculator


464


. The multiplier


462


and the information calculator


464


are used in online (normal) mode. The multiplier


466


and the information calculator


468


are provided for off-line checking.




The signal ε


a


from the output of the adder


490


is provided to a first input of the multiplier


462


and to a second input of the multiplier


462


. An output of the multiplier


462


, being a signal ε


a




2


, is provided to an input of the information calculator


464


. The information calculator


464


computes H


a


(y)≦I(x


a


,y


a


). An output of the information calculator


464


is an information criteria for accuracy and reliability, I(x


a


,y


a


)→max.




The signal ε


b


from the output of the adder


491


is provided to a first input of the multiplier


466


and to a second input of the multiplier


466


. An output of the multiplier


466


, being a signal ε


b




2


, is provided to an input of the information calculator


468


. The information calculator


468


computes H


b


(y)≦I(x


b


,y


b


). An output of the information calculator


468


is an information criteria for accuracy and reliability, I(x


b


,y


b


)→max.




The optimizer


440


includes a second Genetic Algorithm (GA


2


)


444


and an entropy model


442


. The signal I(x


a


,y


a


)→max from the information calculator


464


is provided to a first input of the algorithm (GA


2


)


444


in the optimizer


440


. An entropy signal S→min is provided from an output of the thermodynamic model


442


to a second input of the algorithm GA


2




444


. The signal I(x


b


,y


b


)→max from the information calculator


468


is provided to a third input of the algorithm (GA


2


)


444


in the optimizer


440


.




The signals I(x


a


,y


a


)→max and I(x


b


,y


b


)→max provided to the first and third inputs of the algorithm (GA


2


)


444


are information criteria, and the entropy signal S(k


2


)→min provided to the second input of the algorithm GA


2




444


is a physical criteria based on entropy. An output of the algorithm GA


2




444


is a teaching signal for a FNN


2




486


described below.




The reduced control system


480


includes a reduced sensor set


482


, a control object model


484


, the FNN


2




486


, and a control object


488


. When run in a special off-line checking (verification) mode, the sensor system


482


also includes the set of sensors k. The control object model


484


and the FNN


2




486


together comprise a reduced fuzzy control unit


485


. An output of the control object


488


is provided to an input of the sensor set


482


. An I


2


output of the sensor set


482


contains information from the set n of sensors, where n=(k


1


+k


2


)<m, such that I


2


<I


1


. The information I


2


is provided to a tuning input of the FNN


2




486


, to an input of the control object model


484


, and to an input of the entropy model


442


. The teaching signal


443


from the algorithm GA


2




444


is provided to a teaching signal input of the FNN


2




486


. A control output from the FNN


2




486


is the signal y


a


, which is provided to a control input of the control object


488


.




The control object models


424


and


484


may be either full or partial models. Although

FIGS. 4A and 4B

show the optimal system


420


and the reduced system


480


as separate systems, typically the systems


420


and


480


are the same system. The system


480


is “created” from the system


420


by removing the extra sensors and training the neural network. Thus, typically, the control object models


424


and


484


are the same. The control object


428


and


488


are also typically the same.





FIG. 4B

shows arrow an off-line tuning arrow mark from GA


1




427


to the FNN


1




426


and from the GA


2




444


to the FNN


2




486


.

FIG. 4B

also shows an online learning arrow mark


429


from the sensor system


422


to the FNN


1




426


. The Tuning of the GA


2




444


means to change a set of bonding coefficients in the FNN


2




486


. The bonding coefficients are changed (using, for example, an iterative or trial and error process) so that I(x,y) trends toward a maximum and S trends toward a minimum. In other words, the information of the coded set of the coefficients is sent to the FNN


2




486


from the GA


2




444


as I(x,y) and S are evaluated. Typically, the bonding coefficients in the FNN


2




486


are tuned in off-line mode at a factory or service center.




The teaching signal


429


is a signal that operates on the FNN


1




426


during operation of the optimal control system


420


with an optimal control set. Typically, the teaching signal


429


is provided by sensors that are not used with the reduced control system


480


when the reduced control system is operating in online mode. The GA


1




427


tunes the FNN


1




426


during off-line mode. The signal lines associated with x


b


and y


b


are dashed to indicate that the x


b


and y


b


signals are typically used only during a special off-line checking (i.e., verification) mode. During the verification mode, the reduced control system is run with an optimal set of sensors. The additional sensor information is provided to the optimizer


440


and the optimizer


440


verifies that the reduced control system


480


operates with the desired, almost optimal, accuracy.




For stable and unstable control objects with a non-linear dissipative mathematical model description and reduced number of sensors (or a different set of sensors) the control system design is connected with calculation of an output accuracy of the control object and reliability of the control system according to the information criteria I(x


a


,y


a


)→max and I(x


b


,y


b


)→max. Control system design is also connected with the checking of stability and robustness of the control system and the control object according to the physical criteria S(k


2


)→min.




In a first step, the genetic algorithm GA


2




444


with a fitness function described as the maximum of mutual information between and optimal control signals x


a


and a reduced control signal y


a


is used to develop the teaching signal


443


for the fuzzy neural network FNN


2




486


in an off-line simulation. The fuzzy neural network FNN


2




486


is realized using the learning process with back-propagation for adaptation to the teaching signal, and to develop a lookup-table for changing the parameters of a PID-controller in the controller


485


. This provides sufficient conditions for achieving the desired control reliability with sufficient accuracy.




In a second step, the genetic algorithm GA


2




444


, with a fitness function described as the minimum entropy, S, production rate dS/dt, (calculated according to a mathematical model of the control object


488


), is used to realize a node correction lookup-table in the FNN


2




486


. This approach provides stability and robustness of the reduced control system


480


with reliable and sufficient accuracy of control. This provides a sufficient condition of design for the robust intelligent control system with the reduced number of sensors.




The first and second steps above need not be done in the listed order or sequentially. In the simulation of unstable objects, both steps are preferably done in parallel using the fitness function as the sum of the physical and the information criteria.




After the simulation of a lookup-table for the FNN


2




486


, the mathematical model of control object


484


is changed on the sensor system


482


for checking the qualitative characteristics between the reduced control system with the reduced number of sensors and the reduced control system with a full complement (optimum number) of sensors. Parallel optimization in the GA


2




444


with two fitness functions is used to realize the global correction of a lookup-table in the FNN


2




486


.




The entropy model


442


extracts data from the sensor information I


2


, to help determine the desired number of sensors for measurement and control of the control object


488


.





FIGS. 4A and 4B

show the general case when the reduction of sensors excludes the measurement sensors in the output of the control object and the calculation comparison of control signal on information criteria is possible. The sensor compensator


460


computes the information criteria as a maximum of mutual information between the two control signals x


a


and y


a


(used as the first fitness function for the GA


2




444


). The entropy model


442


present the physical criteria as a minimum production entropy (used as the second fitness function in the GA


2




444


) using the information from sensors


482


. The output of the GA


2




444


is the teaching signal


447


for the FNN


2




486


that is used on-line to develop the reduced control signal y


a


such that the quality of the reduced control signal y


b


is similar to the quality of the optimal control signal x


a


. Thus, the optimizer


440


provides stability and robustness of control (using the physical criteria), and reliability with sufficient accuracy (using the information criteria).




With off-line checking, the optimizer


440


provides correction of the control signal y


a


from the FNN


2




486


using new information criteria. Since the information is additive, it is possible to do the online/off-line steps in sequence, or in parallel. In off-line checking, the sensor system


482


typically uses all of the sensors only for checking of the quality and correction of the control signal y


a


. This approach provides the desired reliability and quality of control, even when the control object


488


is unstable.




For the reduced sensor system


482


(with n sensors) the FNN


2




486


preferably uses a learning and adaptation process instead of a Fuzzy Controller (FC) algorithm.

FIG. 5

illustrates the similarity between a FC


501


and a FNN


540


.




As shown in

FIG. 5A

, the structure of the FNN


540


is similar to the structure of the FC


501


.

FIG. 5A

shows a PID controller


503


having a first input


535


to receive a desired sensor output value and a second input


536


to receive an actual sensor output value


536


. The desired sensor output value is provided to a plus input of an adder


520


and the actual sensor output value


536


is provided to a minus input of the adder


520


. An output of the adder


520


is provided to a proportional output


531


, to an input of an integrator


521


, and to an input of a differentiator


522


. An output of the integrator is provided to an integral output


532


, and an output of the differentiator is provided to a differentiated output


533


.




The Fuzzy logic Controller (FC)


501


includes a fuzzification interface


504


(shown in FIG.


5


B), a Knowledge Base (KB)


502


, a fuzzy logic rules calculator


506


, and a defuzzification interface


508


. The PID outputs


531


-


533


are provided to PID inputs of the fuzzification interface


504


. An output of the fuzzification interface


504


is provided to a first input of the fuzzy logic rules calculator


506


. An output of the KB


502


is provided to a second input of the fuzzy logic rules calculator


506


. An output of the fuzzy logic rules calculator


506


is provided to an input of the defuzzification interface


508


, and an output of the defuzzification interface


508


is provided as a control output U


FC


of the FC


501


.




The Fuzzy Neural Network includes six neuron layers


550


-


555


. The neuron layer


550


is an input layer and has three neurons, corresponding to the three PID output signals


531


-


533


. The neuron layer


550


is an output layer and includes a single neuron. The neuron layers


551


-


553


are hidden layers (e.g., not visible from either the input or the output). The neuron layer


551


has six neurons in two groups of three, where each group is associated with one of the three input neurons in the neuron layer


550


. The neuron layer


552


includes eight neurons. Each neuron in the neuron layer


552


has three inputs, where each input is associated with one of the three groups from the neuron layer


551


. The neuron layer


553


has eight neurons. Each neuron in the neuron layer


553


receives input from all of the neurons in the neuron layer


552


. The neuron layer


554


has is a visible layer having 8 neurons. Each neuron in the neuron layer


554


has four inputs and receives data from one of the neurons in the neuron layer


553


(i.e., the first neuron in the neuron layer


554


receives input from the first neuron in the neuron layer


553


, the second from the second, etc.). Each neuron in the neuron layer


554


also receives input from all of the PID input signals


531


-


533


. The output neuron layer


555


is a single neuron with eight inputs and the single neuron receives data from all of the neurons in the layer


554


(a generalized Takagi-Sugeno FC model).




The number of hidden layers in the FNN


540


typically corresponds to the number of elements in the FC


501


. One distinction between the FC


501


and the FNN


540


lies in the way in which the knowledge base


502


is formed. In the FC


501


, the parameters and membership functions are typically fixed and do not change during of control. By contrast, the FNN


540


uses the learning and adaptation processes, such that when running on-line (e.g., normal mode) the FNN


540


changes the parameters of the membership function in the “IF . . . THEN” portions of the production rules.




The entropy-based control system in

FIGS. 4A and 4B

achieves a desired control accuracy in the output signal of a reduced control system y, as related to an optimal system x, by using the statistical criteria of the accuracy given above as:






M[(x−y)


2


]≦ε


2


  (23)






The optimal output signal x has a constant dispersion σ


x




2


and the variable output signal y has a variable dispersion σ


y




2


. For a desired control accuracy ε, the ε-entropy H


ε


(y) is a measure of the uncertainty in the measurement of the process y as related to the process x, where the processes x and y differ from each other in some metric by ε.




The ε-entropy H


ε


(y) can be calculated as











H
ε



(
y
)


=


H


(
y
)


-

ε





log







n
-
1

ε


-


(

1
-
ε

)






log






1

1
-
ε








(
24
)













where H(y) is the entropy of the process y and n is a number of point measurements in the process y. Further,










H


(
y
)


=

-




k
=
1

n








p
k


log






p
k








(
25
)













where p


k


is the probability that one of the output values of the process y is equal to the value of the output values of the process x.




In general, the form of the ε-entropy is defined as











H
ε



(
y
)


=


inf

p


(

x
,
y

)





I
(

x
,
y

)






(
26
)













with the accuracy M[(x−y)


2


]≦ε


2


.




In general, the amount of information, I(x,y), is given by










I


(

x
,
y

)


=

M


[

log



p


(

x
,
y

)




p


(
x
)




p


(
y
)





]






(
27
)













When x and y are Gaussian process, then










I


(

x
,
y

)


=



1
2



log


(

1
+

1


(

2
+
σ

)


σ



)







where





σ

=


σ
y
2


σ
x
2







(
28
)













The amount of information, given by I(x,y), satisfies the in equality




 I(x,y)≧H


ε


(y)  (29)




and, it may be seen that






M[(x−y)


2


]≧[ε


y


(I(x,y))]


2


  (30)






where the function ε


y


(I(x,y)) is the inverse of H


ε


(y)






H


ε






y






(H)


(y)=H  (31)






and thus






ε


y


(I(x,y))≦ε


y


(H


ε


(y))=M[(x−y)


2


]  (32)






Thus, it is possible to bound the mean-square error M[(x−y)


2


]≦


2


from below by the known amount of information I(x,y).




The calculation of










max

σ
y
2




I


(

x
,
y

)






(
33
)













determines the reliability of the measurements of the process y, as related to the optimal process x, with the accuracy






ε


y


(I(x,y))≦M[(x−y)


2


]≦ε


2


  (34)






In equation (28), the amount of information I(x,y) in the process y as related to the process x will have a maximal value I


max


when σ


y




2


→σ


y




2


. Moreover, at I


max


, the probability p(x,y) of the similarity between the processes x and y is (p(y≈x)→1), and thus, the minimum value of ε


2


is observed. The calculation of equation (33) for the control process y compensates for the loss of information, and minimized the statistical risk, caused by the reduction in the number of sensors.




Engine Control




In one embodiment, the reduced sensor, entropy-based, control system is applied to an engine control system, such as, for example, an internal combustion engine, a piston engine, a diesel engine, a jet engine, a gas-turbine engine, a rocket engine, etc.

FIG. 6

shows an internal combustion piston engine having four sensors, an intake air temperature sensor


602


, a water temperature sensor


604


, a crank angle sensor


606


, and an air-fuel ratio sensor


608


. The air temperature sensor


602


measures the air temperature in an intake manifold


620


. A fuel injector


628


provides fuel to the air in the intake manifold


620


. The intake manifold


620


provides air and fuel to a combustion chamber


622


. Burning of the fuel-air mixture in the combustion chamber


622


drives a piston


628


. The piston


628


is connected to a crank


626


such that movement of the piston


628


turns the crank


626


. The crank angle sensor


606


measures the rotational position of the crank


626


. The water temperature sensor measures the temperature of water in a water jacket


630


surrounding the combustion chamber


622


and the piston


628


. Exhaust gasses from the combustion chamber


622


are provided to an exhaust manifold


624


and the air-fuel ratio sensor


608


measures the ratio of air to fuel in the exhaust gases.





FIG. 7

is a block diagram showing a reduced control system


780


and an optimal control system


720


. The optimal control system


720


, together with an optimizer


740


and a sensor compensator


760


are used to teach the reduced control system


780


. In

FIG. 7

, a desired signal (representing the desired engine output) is provided to an input of the optimal control system


720


and to an input of the reduced control system


780


. The optimal control system


720


, having five sensors, provides an optimal control signal x


a


and an sensor output signal x


b


. A reduced control system


780


provides a reduced control output signal y


a


, and an output sensor signal y


b


. The signals x


b


and y


b


include data from the A/F sensor


608


. The signals x


b


and y


b


are provided to first and second inputs of a subtractor


791


. An output of the subtractor


791


is a signal ε


b


where ε


b


=x


b


−y


b


. The signal ε


b


is provided to a sensor input of a sensor compensator


760


. The signals x


a


and y


a


are provided to first and second inputs of a subtractor


790


. An output of the subtractor


790


is a signal ε


a


where ε=x


a


−y


a


. The signal ε


a


is provided to a control signal input of the sensor compensator


760


. A control information output of the sensor compensator


760


is provided to a control information input of the optimizer


740


. A sensor information output of the sensor compensator


760


is provided to a sensor information input of the optimizer


740


. A sensor signal


783


from the reduced control system


780


is also provided to an input of the optimizer


740


. An output of the optimizer


740


provides a teaching signal


747


to an input of the reduced control system


780


.




The output signal x


b


is provided by an output of a sensor system


722


having five sensors, including, the intake air temperature sensor


602


, the water temperature sensor


604


, the crank angle sensor


606


, the pressure sensor


607


, and the air-fuel ratio (A/F) sensor


608


. The information from the sensor system


722


is a group of signals having optimal information content I


1


. In other words, the information I


1


is the information from the complete set of five sensors in the sensor system


722


.




The output signal x


a


is provided by an output of a control unit


725


. The control signal x


a


is provided to an input of an engine


728


. An output of the engine


728


is provided to an input of the sensor system


722


. Information I


k


from the A/F sensor


608


is provided to an online learning input of a fuzzy neural network (FNN)


726


and to an input of a first Genetic Algorithm (GA


1


)


727


. Information I


k1


from the set of four sensors excluding the A/F sensor


608


is provided to an input of an engine model


724


. An off-line tuning signal output from the algorithm GA


1




727


is provided to an off-line tuning signal input of the FNN


726


. A control output from the FNN


726


is a fuel injector the control signal U


1


, which is provided to a control input of the engine


728


. The signal U


1


is also the signal x


a


. The engine model


724


and the FNN


726


together comprise an optimal control unit


725


.




The sensor compensator


760


includes a multiplier


762


, a multiplier


766


, an information calculator


768


, and an information calculator


764


. The multiplier


762


and the information calculator


764


are used in online (normal) mode. The multiplier


766


and the information calculator


768


are provided for off-line checking.




The signal ε


a


from the output of the adder


790


is provided to a first input of the multiplier


762


and to a second input of the multiplier


762


. An output of the multiplier


762


, being a signal ε


a




2


, is provided to an input of the information calculator


764


. The information calculator


764


computes H


a


(y)≦I(x


a


,y


a


). An output of the information calculator


764


is an information criteria for accuracy and reliability, I(x


a


,y


a


)→max.




The signal ε


b


from the output of the adder


791


is provided to a first input of the multiplier


766


and to a second input of the multiplier


766


. An output of the multiplier


764


, being a signal ε


b




2


, is provided to an input of the information calculator


768


. The information calculator


768


computes H


b


(y)≦I(x


b


,y


b


). An output of the information calculator


768


is an information criteria for accuracy and reliability, I(x


b


,y


b


)→max.




The optimizer


740


includes a second Genetic Algorithm (GA


2


)


744


and a thermodynamic (entropy) model


742


. The signal I(x


a


,y


a


)→max from the information calculator


764


is provided to a first input of the algorithm (GA


2


)


744


in the optimizer


740


. An entropy signal S→min is provided from an output of the thermodynamic model


742


to a second input of the algorithm GA


2




744


. The signal I(x


b


,y


b


)→max from the information calculator


768


is provided to a third input of the algorithm (GA


2


)


744


in the optimizer


740


.




The signals I(x


a


,y


a


)→max and I(x


b


,y


b


)→max provided to the first and third inputs of the algorithm (GA


2


)


744


are information criteria, and the entropy signal S(k


2


)→min provided to the second input of the algorithm GA


2




744


is a physical criteria based on entropy. An output of the algorithm GA


2




744


is a teaching signal for the FNN


786


.




The reduced control system


780


includes a reduced sensor system


782


, an engine model


744


, the FNN


786


, and an engine


788


. The reduced sensor system


782


includes all of the engine sensors in the sensor system


722


except the A/F sensor


608


. When run in a special off-line checking mode, the sensor system


782


also includes the A/F sensor


608


. The engine model


744


and the FNN


786


together comprise a reduced control unit


785


. An output of the engine


788


is provided to an input of the sensor set


782


. An I


2


output of the sensor set


782


contains information from four engine sensors, such that I


2


<I


1


. The information I


2


is provided to an input of the control object model


784


, and to an input of the thermodynamic model


742


. The teaching signal


747


from the algorithm GA


2




744


is provided to a teaching signal input of the FNN


786


. A control output from the FNN


786


is an injector control signal U


2


, which is also the signal y


a


.




Operation of the system shown in

FIG. 7

is in many respects similar to the operation of the system shown in

FIGS. 4A

,


4


B, and


5


B.




The thermodynamic model


742


is built using the thermodynamic relationship between entropy production and temperature information from the water temperature sensor


604


(T


W


) and the air temperature sensor


602


(T


A


). The entropy production S(T


W


,T


A


), is calculated using the relationship









S
=

c
[

ln


(


T
W


T
A


)





]
2


Δτ
-

ln


(


T
W


T
A


)









(
35
)













where Δτ is the duration of a finite process.




The engine model


744


is analyzed from a thermodynamic point of view as a steady-state model where the maximum work, which is equal to the minimum entropy production, is delivered from a finite resource fluid and a bath. Moreover, the engine model is analyzed as a dissipative, finite-time, generalization of the evolutionary Carnot problem in which the temperature driving force between two interacting subsystems varies with the constant time Δτ (a reversible cycle).




The presence of reversible cycles fixes, automatically, the first-law efficiency of each infinitesimal stage at the Carnot level









η
=

1
-


T
e

T






(
36
)













where T is the instantaneous temperature of a finite resource (such as the water temperature) and T


e


is the temperature of the environment or an infinite bath (such as the air temperature). The unit mass of a resource releases the heat dq=−cdT, where c is the specific heat. The classical rate and duration-dependent function of available energy (dissipative energy) E


x


follows by integration of the product











-
c






η







T


=


-

c


(

1
-


T
e

T


)






T






(
37
)













between the limits T and T


e


. The integration yields the well known expression









W
=




T

T
e





c


(

1
-


T
e

T


)









T



=


Δ





h

-


T
e


Δ





s







(
38
)













The problem becomes non-trivial in the case where a finite-rate process is considered, because in the finite-rate process, the efficiency differs from the Carnot efficiency. The finite-rate efficiency should be determined before the integration of the product −cdT can be evaluated. The integration then leads to a generalized available energy associated with the extremal release of the mechanical work in a finite time.




The local efficiency of an infinitesimal unit of the finite-rate process is









η
=



W




Q
1







(
39
)













where Q


1


is the cumulative heat flux flowing from the upper reservoir. A flux Q


2


is the cumulative heat flux flowing to the lower reservoir. While this local first-law efficiency is still described by the Carnot formula









η
=

1
-


T
2



T
1








(
40
)













were T


1


≦T


1


′≦T


2


′≦T


2


, this efficiency is nonetheless lower than the efficiency of the unit working between the boundary temperatures T


1


and T


2


=T


e


, as the former applies to the intermediate temperatures T


1


′ and T


2


′. By solving equation (40) along with the reversible entropy balance of the Carnot differential subsystem


















γ
1



(


T
1

-

T
1



)




T
1



=









γ
2



(


T
2

-

T
2



)




T
2







(
41
)













one obtains the primed temperatures as functions of the variables T


1


, T


2


=T


e


and η. In equation (41), the parameters γ


1


and γ


2


are the partial conductances, and thus link the heat sources with the working fluid of the engine at high and low temperatures. The differentials of the parameters γ


1


and γ


2


can be expressed as dγ


1





i


dA


i


, for i=1,2, where the values α


i


are the heat transfer coefficients and the values dA


i


are the upper and lower exchange surface areas.




The associated driving heat flux dQ


1


=dγ


1


(T


1


−T


1


′) is then found in the form
















Q
1


=







γ


[


T
1

-


1

(

1
-
η

)




T
2



]







(
42
)













from which the efficiency-power characteristic follows in the form









η
=

1
-


T
2



T
1

-




Q
1








γ









(
43
)













From the energy balance of the driving fluid, the heat power variable Q


1


satisfies dQ=−GcdT, where dT is the differential temperature drop of the driving fluid, G is the finite-mass flux of the driving fluid whose constant specific heat capacity equals c. Using the above definitions, and the differential heat balance of the driving fluid, the control term dQ


1


/dγ from equation (43) may be written (with the subscript 1 omitted) as












Q



γ


=



T



τ






(
44
)













The negative of the derivative dQ/dγ is the control variable u of the process. In short, the above equation says that u={dot over (T)}, or that the control variable u equals the rate of the temperature change with respect to the non-dimensional time τ. Thus, equation (43) becomes the simple, finite-rate generalization of the Carnot formula









η
=

1
-


T
e


T
+

T
.








(
45
)













When T>T


e


then the derivative {dot over (T)} is negative for the engine mode. This is because the driving fluid releases the energy to the engine as part of the work production.




Work Functionals for an Infinite Sequence of Infinitesimal Processes




The cumulative power delivered per unit fluid flow W/G is obtained by the integration of the product of n and dQ/G=−cdT between an arbitrary initial temperature T


i


and a final temperature T


f


of the fluid. This integration yields the specific work of the flowing fluid in the form of the functional











W

[


T
i

,

T
f


]




W
G


=

-




T
i


T
f





c


(

1
-


T
e


T
+

T
.




)




T
.








τ








(
46
)













The notation [T


i


,T


f


] means the passage of the vector T≡(T,τ) from its initial state to its final state. For the above functional, the work maximization problem can be stated for the engine mode of the process by













W

m





a





x


=





max


{

-




T
i


T
f





L


(

T
,

T
.


)









τ




}








=





max


{

-




T
i


T
f





c


(

1
-


T
e


T
+

T
.




)




T
.








τ




}









(
47
)













In equation (47), the function L(T,{dot over (T)}) is the Lagrangian. The above Lagrangian functional represents the total power per unit mass flux of the fluid, which is the quantity of the specific work dimension, and thus the direct relation to the specific energy of the fluid flow. In the quasi-static limit of vanishing rates, where dT/dτ=0, the above work functional represents the change of the classical energy according to equation (38).




For the engine mode of the process, the dissipative energy itself is obtained as the maximum of the functional in equation (47), with the integration limits T


i


=T


e


and T


f


=T.




An alternative form of the specific work can be written as the functional














W

[


T
i

,

T
f


]




W
G


=









T
i


T
f





-

c


(

1
-


T
e


T
.



)






τ



-


T
e






T
i


T
f




c




T
.

2


T


(

T
+

T
.


)













=






-




T
i


T
f





c


(

1
-


T
e


T
.



)





T




-


T
e


S









(
48
)













in which the first term is the classical “reversible” term and the second term is the product of the equilibrium temperature and the entropy production S, where









S
=




T
i


T
e




c




T
.

2


T


(

T
+

T
.


)










τ







(
49
)













The entropy generation rate is referred here to the unit mass flux of the driving fluid, and thus to the specific entropy dimension of the quantity S. The quantity S is distinguished from the specific entropy of the driving fluid, s in equation (38). The entropy generation S is the quadratic function of the process rate, u=dT/dτ, only in the case when the work is not produced, corresponding to the vanishing efficiency η=0 or the equality T+{dot over (T)}=T


e


. For the activity heat transfer (non-vanishing η) case, the entropy production appears to be a non-quadratic function of rates, represented by the integrand of equation (48). Applying the maximum operation for the functional in equation (48) at the fixed temperatures and times, it is seen that the role of the first term (i.e., the potential term) is inessential, and the problem of the maximum released work is equivalent to the associated problem of the minimum entropy production. This confirms the role of the entropy generation minimization in the context of the extremum work problem. The consequence of the conclusion is that a problem of the extremal work and an associated fixed-end problem of the minimum entropy generation have the same solutions.




Hamilton-Jacobi Approach to Minimum Entropy Generation




The work extremization problems can be broken down to the variational calculus for the Lagrangian









L
=


c


(

1
-


T
e



T
.

-
T



)




T
.






(
50
)













The Euler-Lagrangian equations for the problem of extremal work and the minimum entropy production lead to a second-order differential equation






T{umlaut over (T)}−{dot over (T)}


c


=0  (51)






which characterizes the optimal trajectories of all considered processes. A first integral of the above equation provides a formal analog of the mechanical energy as










E





L




T
.



-
L


=

c








T
e




T
.

2




(


T
.

+
T

)

2







(
52
)













An equation for the optimal temperature flow from the condition E=h in equation (38) is given by










T
.

=





±
T




h

cT
e





1
-

±


h

cT
e








ξ





T





where





ξ


=


ln



T
f


T
i





τ
f

-

τ
i








(
53
)













The power of dynamic programming (DP) methods when applied to problems of this sort lies in the fact that regardless of local constraints, functions satisfy an equation of the Hamilton-Jacobi-Bellman (HJB-equation) with the same state variables as those for the unconstrained problem. Only the numerical values of the optimizing control sets and the optimal performance functions differ in constrained and unconstrained cases.




The engine problem may be correctly described by a backward HJB-equation. The backward HJB-equation is associated with the optimal work or energy as an optimal integral function I defined on the initial states (i.e., temperatures), and, accordingly, refers to the engine mode or process approaching an equilibrium. The function I is often called the optimal performance function for the work integral.




The maximum work delivery (whether constrained or unconstrained) is governed by the characteristic function













I


(


τ
f

,

T
f

,

τ
i

,

T
i


)








max






W

[


T
i

,

T
f


]









=





max


{




T
i


T
f





[


-

c


(

1
-


T
e


T
+
u



)




u

]








t



}









(
54
)













The quantity I in equation (54) describes the extremal value of the work integral W(T


i


,T


f


). It characterizes the extremal value of the work released for the prescribed temperatures T


i


and T


f


when the total process duration is τ


ƒ


−τ


i


. While the knowledge of the characteristic function I only is sufficient for a description of the extremal properties of the problem, other functions of this sort are nonetheless useful for characterization of the problem.




Here the problem is transformed into the equivalent problem in which one seeks the maximum of the final work coordinate x


o




ƒ


=W


ƒ


for the system described by the following set of differential equations















W



τ


=



-

c


(

1
-


T
e


u
+
T



)




u




f
0



(

T
,
u

)












T



τ


=

u



f
1



(

T
,
u

)













x
2




τ


=

1



f
2



(

T
,
u

)










(
55
)













The state of the system described by equations (55) is described by the enlarged state vector X=(x


0


, x


1


, x


2


) which is composed of the three state coordinates x


0


=W, x


1


=T, and x


2


=τ. At this point, it is convenient to introduce a set of optimal performance functions Θ


i


and Θ


f


corresponding to the initial and final work coordinates, respectively. The function Θ


i


works in a space of one dimension larger than I, involves the work coordinate x


0


=W, and is defined as












max
u



W
f





Θ
i



(


W
i

,

τ
i

,

T
i

,

τ
f

,

T
f


)



=


W
i

+

I


(


τ
i

,

T
i

,

τ
f

,

T
f


)







(
56
)













A function V is also conveniently defined in the enlarged space of variables as






V=W


ƒ


−Θ


i


(W


i





i


,T


i





ƒ


, T


ƒ


)=W


ƒ


−W


i


−I(τ


i


,T


i





ƒ


,T


ƒ


)  (57)






Two mutually equal maxima at the constant W


i


and at the constant W


f


are described by the extremal functions






V


i


(W


i





i


,T


i





ƒ


,T


ƒ


)=V


ƒ





i


,T


i


,W


ƒ





ƒ


,T


ƒ


)=0  (58))






which vanish identically along all optimal paths. These functions can be written in terms of a wave-function V as follows:




 max V=max{W


ƒ


−W


i


−I(τ


ƒ


,T


ƒ





i


,T


i


)}=0  (59)




The equation for the integral work function I in the narrow space of the coordinates (T,τ), which does not involve the coordinates W, follows immediately from the condition V=0. (Development of the backward DP algorithm and equations for Θ


i


(X


i


) are provided later in the discussion that follows.)




Solution of the HJB equations begins by noting that the system of three state equations in equation (55) have the state variables x


0


=W, x


1


=T and x


2


=τ. The state equations from equation (55) can be in a general form as














x
β




τ


=


f
β



(

x
,
u

)



,

β
=
0

,
1
,
2




(
60
)













For small Δτ, from equation (60) it follows that






Δx


β





β


(x,u)Δτ+O(ε


2


)  (61)






In equation (61), the symbol O(ε


2


) means second-order and higher terms. The second-order and higher terms posses the property that











lim


Δ





τ


0





O


(

ε
2

)



Δ





τ




0.




(
62
)













The optimal final work for the whole path in the interval [τ


i





f


] is the maximum of the criterion






W


ƒ





i


(x


i


+Δx)=Θ


i


(W


i


+ΔW,T


i


+ΔT,τ


i


+Δτ)  (63)






Expanding equation (63) in a Taylor series yields













W
f

=







Θ
i



(

X
i

)


+





Θ
i





X
β
i




Δ






X
β


+

O


(

ε
2

)









=







Θ
i



(


W
i

,

T
i

,

τ
i


)


+





Θ
i





W
i




Δ





W

+





Θ
i





T
i




Δ





T





+





Θ
i





τ
i




Δ





τ

+

O


(

ε
2

)










(
64
)













Substituting equation (62) into equation (64) and performing an appropriate extremization in accordance with Bellman's principle of optimality, yields, for the variation of the initial point











max

u
i




W
f


=


max

u
i




{



Θ
i



(

X
i

)


+





Θ
i





X
β
i






f
β



(

X
,
u

)



Δ





τ

+

O


(

ε
2

)



}






(
65
)













After reduction of Θ


i


, and the division of both sides of equation (65) by Δτ, the passage to the limit Δτ→0 subject to the condition











lim


Δ





τ


0





O


(

ε
2

)



Δ





τ




0




(
66
)













yields the backward HJB equation of the optimal control problem. For the initial point of the extremal path, the backward DP equation is










max



{





Θ
i





X
β
i






f
β



(

X
,
u

)



}


u
i



=


max



{






Θ
i





W
i







W
.

i



(


T
i

,

u
i


)



+





Θ
i





T
i







T
.

i



(


T
i

,

u
i


)



+




Θ
i





τ
i




}




u
i




=


max



{




Θ
i






τ
i




}


u
i



=


-


min

u
i




{



V




τ
i



}



=


max


{



V




(

-

τ
i


)



}


=
0








(
67
)













The properties of V=W


i


−Θ


i


have been used in the second line of the above equation.




The partial derivative of V with respect to the independent variable τ can remain outside of the bracket of this equation as well. Using the relationships












V




W
i



=


-




Θ
i





W
i




=

-
1






(
68
)













and {dot over (W)}=ƒ


0


=−L (see e.g., equations (47) and (50)), leads to the extremal work function













V




τ
i



+


min

u
i




{





V




T
i





u
i


+


L
i



(


T
i

,

u
i


)



}



=

0






(

max






W
f


)






(
69
)













In terms of the integral function of optimal work, I=W


ƒ


−W


i


−V, equation (69) becomes













I




τ
i



+


min

u
i




{





I




T
i





u
i


+


f
0



(


T
i

,

u
i


)



}



=
0




(
70
)













As long as the optimal control u is found in terms of the state, time, and gradient components of the extremal performance function I, the passage from the quasi-linear HJB equation to the corresponding nonlinear Hamilton-Jacobi equation is possible.




The maximization of equation (70) with respect to the rate it leads to two equations, the first of which describes the optimal control u expressed through the variables T and z=∂I/∂T as follows:












I



T


=

-









f
0



(

T
,
u

)





u







(
71
)













and the second is the original equation (70) without the extremization sign













I



τ


+




I



T



u

+


f
0



(

T
,
u

)



=
0




(
72
)













In the proceeding two equations, the index i is omitted. Using the momentum-type variable z=∂I/∂T, and using equation (71) written in the form









z
=


-





f
0



(

T
,
u

)





u



=




L


(

T
,
u

)





u







(
73
)













leads to the energy-type Hamiltonian of the extremal process as






H=zu(z,T)+ƒ


0


(z,T)  (74)






Using the energy-type Hamiltonian and equation and equation (72) leads to the Hamilton-Jacobi equation for the integral I as













I



τ


+

H


(

T
,



I



τ



)



=
0




(
75
)













Equation (75) differs from the HJB equation in that equation (75) refers to an extremal path only, and H is the extremal Hamiltonian. The above formulas are applied to the concrete Lagrangian L=−ƒ


0


, where ƒ


0


is the intensity of the mechanical work production.




The basic integral W


[T






i






,T






ƒ






]


written in the form of equation (47)













W

m





a





x


=





max


{

-




T
i


T
f





L


(

T
,

T
.


)









τ




}








=





max


{

-




τ
i


τ
f





c


(

1
-


T
e


T
+

T
.




)




T
.








τ




}









(
76
)













whose extremal value is the function I(T


i





i


,T


ƒ





ƒ


). The momentum-like variable (equal to the temperature adjoint) is then










z





f
0




u



=

c


(

1
-



T
e


T



(

T
+
u

)

2



)






(
77
)





and










u
=



(



T
e


T


1
-

z
c



)


1
/
2


-
T





(
78
)













The Hamilton-Jacobi partial differential equation for the maximum work problem (engine mode of the system) deals with the initial coordinates, and has the form













I



τ


+

H


(

T
,



I



τ



)



=
0




(
79
)







H


(

T
,



I



T



)


=


c


[



T
e


-


T


(

1
-


1
c





I



T




)




]


2



















The variational front-end problem for the maximum work W is equivalent to the variational fixed-end problem of the minimum entropy production.




Hamilton-Jacobi Approach to Minimum Entropy Production




The specific entropy production is described by the functional from equation (49) as










S
σ

=




0

τ
f





L
σ








τ



=



0

τ
f




c



u
2


T


(

T
+
u

)










τ








(
80
)













The minimum of the functional can be described by the optimal function I


σ


(T


i


,mτ


i


,T


ƒ





ƒ


) and the Hamilton-Jacobi equation is then














I
σ




τ


+

H
σ


=
0




(
81
)







H
σ

=


c


[

1
-


1
-


1
c


T





I
σ




T






]


2



















The two functionals (that of the work and that of the entropy generation) yield the same is extremal.




Considering the Hamilton-Jacobi equation for the engine, from equation (54) integration along an extremal path leads to a function which describes the optimal specific work










I


(


T
i

,

T
f

,

τ
i

,

τ
f


)


=


c


(


T
f

-

T
i


)


-



T
e


1
+
ξ



ln







T
i


T
f








(
82
)













where the parameter ξ is defined in equation (53). The extremal specific work between two arbitrary states follows in the form













I


(


T
i

,

T
f

,

τ
i

,

τ
f


)


=






c


(


T
i

-

T
f


)


-


cT
e


ln



T
i


T
f



-















cT
e



[

ln



T
i


T
f



]


2



τ
f

-

τ
i

-

ln



T
i


T
f












(
83
)













From equation (83), with T


i


=T


w


and T


f


=T


A


, then the minimal integral of the entropy production is











S
σ

=








c


[

ln



T
W


T
A



]


2



Δτ
-

[

ln



T
W


T
A



]




,


where





Δτ

=


τ
f

-

τ
i







(
84
)













The function from equation (83) satisfies the backward Hamilton-Jacobi equation, which is equation (79), and the function (84) satisfies the Hamilton-Jacobi equation (81).




The Lyapunnov function ℑ for the system of equation (60) can be defined as









𝔍
=



1
2






β
=
1

3







x
β
2



+


1
2



S
2







(
85
)





and













𝔍



τ


=





β
=
1

3








x
β



f
β



+


(


S
u

-

S
c


)



(





S
u




τ


-




S
c




τ



)







(
86
)













The relation (86) is the necessary and sufficient conditions for stability and robustness of control with minimum entropy in the control object and the control system.




Suspension Control




In one embodiment, the reduced control system of

FIGS. 4A and 4B

is applied to a suspension control system, such as, for example, in an automobile, truck, tank, motorcycle, etc.





FIG. 8

is a schematic diagram of one half of an automobile suspension system. In

FIG. 8

, a right wheel


802


is connected to a right axle


803


. A spring


804


controls the angle of the axle


803


with respect to a body


801


. A left wheel


812


is connected to a left axle


813


and a spring


814


controls the angle of the axle


813


. A torsion bar


820


controls the angle of the left axle


803


with respect to the right axle


813


.





FIG. 9

is a block diagram of a control system for controlling the automobile suspension system shown in FIG.


8


. The block diagram of

FIG. 9

is similar to the diagram of

FIG. 4B

showing an optimal control system and a reduced control system


980


.




The optimal control system


920


includes an optimal sensor system


922


that provides information I


1


to an input of a control object model


924


. An output of the control object model


924


is provided to a first genetic analyzer GA


1




927


, and the GA


1




927


provides an off-line tuning signal to a first fuzzy neural network FNN


1




926


. An output from the FNN


1




926


, and an output from the optimal sensor system


922


are provided to a control unit


925


. An output from the control unit


925


is provided to a control object


928


. The control object


928


is the automobile suspension and body shown in FIG.


8


.




The reduced control system


980


includes a reduced sensor system


982


that provides reduced sensor information I


2


to an input of a control unit


985


. The control unit


985


is configured by a second fuzzy neural network FNN


2




986


. An output from the control unit


985


is provided to a control object


988


. The control object


988


is also the automobile suspension and body shown in FIG.


8


.




The FNN


2




986


is tuned by a second genetic analyzer GA


2


, that, (like the GA


2




444


in

FIG. 4.

) maximizes an information signal I(x,y) and minimizes an entropy signal S. The information signal I(x,y) is computed from a difference between the control and sensor signals produced by the optimum control system


920


and the reduced control system


980


. The entropy signal S is computed from a mathematical model of the control object


988


.




The optimal sensor system


922


includes a pitch angle sensor, a roll angle sensor, four position sensors, and four angle sensors as described below.




The reduced sensor system


922


includes a single sensor, a vertical accelerometer placed near the center of the vehicle


801


.




A Half Car Coordinate Transformation




1. Description of transformation matrices




1.1 Global reference coordinate x


r


,y


r


,z


y


{r} is assumed to be at the pivot center P


r


of a rear axle.




The following are the transformation matrices to describe the local coordinates for:




{2} the gravity center of the body.




{7} the gravity center of the suspension.




{10} the gravity center of the arm.




{12} the gravity center of the wheel.




{13} the touch point of the wheel to the road.




{14} the stabilizer linkage point.




1.2 Transformation matrices.




Rotating {r} along y


r


with angle β makes a local coordinate system x


0r


, y


0r


, z


0r


{0r} with a transformation matrix


r




0r


T.













0

r

r


T

=

[




cos





β



0



sin





β



0




0


1


0


0






-
sin






β



0



cos





β



0




0


0


0


1



]





(
87
)













Transferring {0r} through the vector (a


1


, 0, 0) makes a local coordinate system x


0f


,y


0f


,z


0f


{0f} with a transformation matrix


0r




0f


T.













0

f


0

r



T

=

[



1


0


0



a
1





0


1


0


0




0


0


1


0




0


0


0


1



]





(
88
)













The above procedure is repeated to create other local coordinate systems with the following transformation matrices.













1

f


0

f



T

=

[



1


0


0


0




0



cos





α





-
sin






α



0




0



sin





α




cos





α



0




0


0


0


1



]





(
89
)









2

1

f



T

=

[



1


0


0



a
0





0


1


0



b
0





0


0


1



c
0





0


0


0


1



]





(
90
)













1.3 Coordinates for the front wheels


802


,


812


(index n: i for the left, ii for the right) are generated as follows. Transferring {1f} through the vector (0, b


2n


, 0) makes local coordinate system x


3n


, y


3n


, z


3n


{3n} with transformation matrix


1f




3n


T.













3

n


1

f



T

=

[



1


0


0


0




0


1


0



b

2

n






0


0


1


0




0


0


0


1



]





(
91
)










4

n


3

n



T

=

[



1


0


0


0




0



cos






γ
n






-
sin







γ
n




0




0



sin






γ
n





cos






γ
n




0




0


0


0


1



]





(
92
)










5

n


4

n



T

=

[



1


0


0


0




0


1


0


0




0


0


1



c

1

n






0


0


0


1



]





(
93
)










6

n


5

n



T

=

[



1


0


0


0




0



cos






η
n






-
sin







η
n




0




0



sin






η
n





cos






η
n




0




0


0


0


1



]





(
94
)










7

n


6

n



T

=

[



1


0


0


0




0


1


0


0




0


0


1



z

6

n






0


0


0


1



]





(
95
)










8

n


4

n



T

=

[



1


0


0


0




0


1


0


0




0


0


1



c

2

n






0


0


0


1



]





(
96
)










9

n


8

n



T

=

[



1


0


0


0




0



cos






θ
n






-
sin







θ
n




0




0



sin






θ
n





cos






θ
n




0




0


0


0


1



]





(
97
)










10

n


9

n



T

=

[



1


0


0


0




0


1


0



e

1

n






0


0


1


0




0


0


0


1



]





(
98
)










11

n


9

n



T

=

[



1


0


0


0




0


1


0



e

3

n






0


0


1


0




0


0


0


1



]





(
99
)










12

n


11

n



T

=

[



1


0


0


0




0



cos






ζ
n






-
sin







ζ
n




0




0



sin






ζ
n





cos






ζ
n




0




0


0


0


1



]





(
100
)










13

n


12

n



T

=

[



1


0


0


0




0


1


0


0




0


0


1



z

12

n






0


0


0


1



]





(
101
)










14

n


9

n



T

=

[



1


0


0


0




0


1


0



e

0

n






0


0


1


0




0


0


0


1



]





(
102
)













1.4 Some matrices are sub-assembled to make the calculation simpler.
















1

f

r


T

=






T

0

f


0

r



0

r

r



T

1

f


0

f



T







=







[




cos





β



0



sin





β



0




0


1


0


0






-
sin






β



0



cos





β



0




0


0


0


1



]



[



1


0


0



a
1





0


1


0


0




0


0


1


0




0


0


0


1



]




[



1


0


0


0




0



cos





α





-
sin






α



0




0



sin





α




cos





α



0




0


0


0


1



]








=






[




cos





β



0



sin





β





a
1


cos





β





0


1


0


0






-
sin






β



0



cos





β





-

a
1



sin





β





0


0


0


1



]



[



1


0


0


0




0



cos





α





-
sin






α



0




0



sin





α




cos





α



0




0


0


0


1



]








=





[




cos





β




sin





βsin





α




sin





β





cos





α





a
1


cos





β





0



cos





α





-
sin






α



0






-
sin






β




cos





β





sin





α




cos





β





cos





α





-

a
1



sin





β





0


0


0


1



]








(
103
)













4

n

r


T

=






T

3

n


1

f



1

f

r



T

4

n


3

n



T







=







[




cos





β




sin





βsin





α




sin





β





cos





α





a
1


cos





β





0



cos





α





-
sin






α



0






-
sin






β




cos





β





sin





α




cos





β





cos





α





-

a
1



sin





β





0


0


0


1



]



[



1


0


0


0




0


1


0



b

2

n






0


0


1


0




0


0


0


1



]




[



1


0


0


0




0



cos






γ
n






-
sin







γ
n




0




0



sin






γ
n





cos






γ
n




0




0


0


0


1



]








=





[




cos





β




sin






βsin


(

α
+

γ
n


)






sin






βcos


(

α
+

γ
n


)








b

2

n



sin





βsinα

+


a
1


cos





β






0



cos


(

α
+

γ
n


)





-

sin


(

α
+

γ
n


)







b

2

n



cos





α







-
sin






β




cos






βsin


(

α
+

γ
n


)






cos






βcos


(

α
+

γ
n


)








b

2

n



cos





βsin





α

-


a
1


sin





β






0


0


0


1



]








(
104
)













7

n


4

n



T

=






T

6

n


5

n



5

n


4

n




T

7

n


6

n



T







=







[



1


0


0


0




0


1


0


0




0


0


1



c

1

n






0


0


0


1



]



[



1


0


0


0




0



cos






η
n






-
sin







η
n




0




0



sin






η
n





cos






η
n




0




0


0


0


1



]




[



1


0


0


0




0


1


0


0




0


0


1



z

6

n






0


0


0


1



]








=






[



1


0


0


0




0



cos






η
n






-
sin







η
n




0




0



sin






η
n





cos






η
n





c

1

n






0


0


0


1



]



[



1


0


0


0




0


1


0


0




0


0


1



z

6

n






0


0


0


1



]








=





[



1


0


0


0




0



cos






η
n






-
sin







η
n






-

z

6

n




sin






η
n






0



sin






η
n





cos






η
n






c

1

n


+


z

6

n



cos






η
n







0


0


0


1



]








(
105
)













10

n


4

n



T

=






T

9

n


8

n



8

n


4

n




T

10

n


9

n



T







=







[



1


0


0


0




0


1


0


0




0


0


1



c

2

n






0


0


0


1



]



[



1


0


0


0




0



cos






θ
n






-
sin







θ
n




0




0



sin






θ
n





cos






θ
n




0




0


0


0


1



]




[



1


0


0


0




0


1


0



e

1

n






0


0


1


0




0


0


0


1



]








=






[



1


0


0


0




0



cos






θ
n






-
sin







θ
n




0




0



sin






θ
n





cos






θ
n





c

2

n






0


0


0


1



]



[



1


0


0


0




0


1


0



e

1

n






0


0


1


0




0


0


0


1



]








=





[



1


0


0


0




0



cos






θ
n






-
sin







θ
n






e

1

n



cos






θ
n






0



sin






θ
n





cos






θ
n






c

2

n


+


e

1

n



sin






θ
n







0


0


0


1



]








(
106
)













12

n


4

n



T

=






T

9

n


8

n



8

n


4

n




T

11

n


9

n




T

12

n


11

n



T







=








[



1


0


0


0




0


1


0


0




0


0


1



c

2

n






0


0


0


1



]



[



1


0


0


0




0



cos






θ
n






-
sin







θ
n




0




0



sin






θ
n





cos






θ
n




0




0


0


0


1



]




[



1


0


0


0




0


1


0



e

3

n






0


0


1


0




0


0


0


1



]




[



1


0


0


0




0



cos






ζ
n






-
sin







ζ
n




0




0



sin






ζ
n





cos






ζ
n




0




0


0


0


1



]








=







[



1


0


0


0




0



cos






θ
n






-
sin







θ
n




0




0



sin






θ
n





cos






θ
n





c

2

n






0


0


0


1



]



[



1


0


0


0




0


1


0



e

3

n






0


0


1


0




0


0


0


1



]




[



1


0


0


0




0



cos






ζ
n






-
sin







ζ
n




0




0



sin






ζ
n





cos






ζ
n




0




0


0


0


1



]








=






[



1


0


0


0




0



cos






θ
n






-
sin







θ
n






e

3

n



cos






θ
n






0



sin






θ
n





cos






θ
n






c

2

n


+


e

3

n



sin






θ
n







0


0


0


1



]



[



1


0


0


0




0



cos






ζ
n






-
sin







ζ
n




0




0



sin






ζ
n





cos






ζ
n




0




0


0


0


1



]








=





[



1


0


0


0




0



cos






(


θ
n

+

ζ
n


)






-
sin







(


θ
n

+

ζ
n


)






e

3

n



cos






θ
n






0



sin






(


θ
n

+

ζ
n


)





cos






(


θ
n

+

ζ
n


)






c

2

n


+


e

3

n



sin






θ
n







0


0


0


1



]








(
107
)













2. Description of all the parts of the model both in local coordinate systems and relations to the coordinate {r} or {1f}.




2.1 Description in local coordinate systems.










P
body
2

=


P

susp
.
n


7

n


=


P

urm
.
n


10

n


=


P

wheel
.
n


12

n


=


P

touchpoint
.
n


13

n


=


P

stab
.
n


14

n


=

[









0




0







0







1



]










(
108
)













2.2 Description in global reference coordinate system {r}.













P
body
r

=






T
2

1

f



1

f

r



TP
body
2








=







[




cos





β




sin





βsin





α




sin





βcos





α





a
1


cos





β





0



cos





α





-
sin






α



0






-
sin






β




cos





β





sin





α




cos





β





cos





α





-

a
1



sin





β





0


0


0


1



]



[



1


0


0



a
0





0


1


0



b
0





0


0


1



c
0





0


0


0


1



]




[









0




0







0







1



]








=





[












a
0


cos





β

+


b
0


sin





βsin





α

+


c
0


sin





βcos





α

+


a
1


cos





β









b
0


cos





α

-


c
0


sin





α












-

a
0



sin





β

+


b
0


cos





βsin





α

+


c
0


cos





βcos





α

-


a
1


sin





β









1



]








(
109
)










P
suspn
r

=






T

7

n


4

n



4

n

r



TP
suspn

7

n









=





[




cos





β




sin






βsin


(

α
+

γ
n


)






sin





β






cos


(

α
+

γ
n


)








b

2

n



sin





β





sin





α

+


a
1


cos





β






0



cos


(

α
+

γ
n


)





-

sin


(

α
+

γ
n


)







b

2

n



cos





α







-
sin






β




cos






βsin


(

α
+

γ
n


)






cos






βcos


(

α
+

γ
n


)








b

2

n



cos





β





sin





α

-


a
1


sin





β






0


0


0


1



]













[



1


0


0


0




0



cos






η
n






-
sin







η
n






-

z

6

n




sin






η
n






0



sin






η
n





cos






η
n






c

1

n


+


z

6

n



cos






η
n







0


0


0


1



]



[









0




0







0







1



]








=





[












{



z

6

n




cos


(

α
+

γ
n

+

n
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


sin





β

+


a
1


cos





β









-

z

6

n





sin


(

α
+

γ
n

+

η
n


)



-


c

1

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α












{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-


a
1


sin





β









1



]








(
110
)










P
armn
r

=






T

10

n


4

n



4

n

r



TP
armn

10

n









=





[




cos





β




sin






βsin


(

α
+

γ
n


)






sin





β






cos


(

α
+

γ
n


)








b

2

n



sin





β





sin





α

+


a
1


cos





β






0



cos


(

α
+

γ
n


)





-

sin


(

α
+

γ
n


)







b

2

n



cos





α







-
sin






β




cos






βsin


(

α
+

γ
n


)






cos






βcos


(

α
+

γ
n


)








b

2

n



cos





β





sin





α

-


a
1


sin





β






0


0


0


1



]













[



1


0


0


0




0



cos


(


θ
n

+

ζ
n


)





-

sin


(


θ
n

+

ζ
n


)







e

3

n



cos






θ
n






0



sin


(


θ
n

+

ζ
n


)





cos


(


θ
n

+

ζ
n


)






c

2

n


+


e

3

n



sin






θ
n







0


0


0


1



]



[









0




0







0







1



]








=





[












{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


sin





β

+


a
1


cos





β









e

1

n




cos


(

α
+

γ
n

+

θ
n


)



-


c

2

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α












{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-


a
1


sin





β









1



]








(
111
)










P

wheel
.
n

r

=






T

12

n


4

n



4

n

r



TP

wheel
.
n


12

n









=





[




cos





β




sin






βsin


(

α
+

γ
n


)






sin





β






cos


(

α
+

γ
n


)








b

2

n



sin





β





sin





α

+


a
1


cos





β






0



cos


(

α
+

γ
n


)





-

sin


(

α
+

γ
n


)







b

2

n



cos





α







-
sin






β




cos






βsin


(

α
+

γ
n


)






cos






βcos


(

α
+

γ
n


)








b

2

n



cos





β





sin





α

-


a
1


sin





β






0


0


0


1



]













[



1


0


0


0




0



cos


(


θ
n

+

ζ
n


)





-

sin


(


θ
n

+

ζ
n


)







e

3

n



cos






θ
n






0



sin


(


θ
n

+

ζ
n


)





cos


(


θ
n

+

ζ
n


)






c

2

n


+


e

3

n



sin






θ
n







0


0


0


1



]



[









0




0







0







1



]








=





[












{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


sin





β

+


a
1


cos





β









e

1

n




cos


(

α
+

γ
n

+

θ
n


)



-


c

2

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α












{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-


a
1


sin





β









1



]








(
112
)










P

touchpoint
.
n

r

=






T

12

n


4

n



4

n

r



T

13

n


12

n




T

touchpoint
.
n


13

n









=





[




cos





β




sin






βsin


(

α
+

γ
n


)






sin






βcos


(

α
+

γ
n


)








b

2

n



sin





β





sin





α

+


a
1


cos





β






0



cos


(

α
+

γ
n


)





-

sin


(

α
+

γ
n


)







b

2

n



cos





α







-
sin






β




cos






βsin


(

α
+

γ
n


)






cos






βcos


(

α
+

γ
n


)








b

2

n



cos





β





sin





α

-


a
1


sin





β






0


0


0


1



]














[



1


0


0


0




0



cos


(


θ
n

+

ζ
n


)





-

sin


(


θ
n

+

ζ
n


)







e

3

n



cos






θ
n






0



sin


(


θ
n

+

ζ
n


)





cos


(


θ
n

+

ζ
n


)






c

2

n


+


e

3

n



sin






θ
n







0


0


0


1



]



[



1


0


0


0




0


1


0


0




0


0


1



z

12

n






0


0


0


1



]




[









0




0







0







1



]








=





[












{



z

12

n



cos





α

+


e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c
n



cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


sin





β

+


a
1


cos





β









-

z

12

n




sin





α

+


e

3

n




cos


(

α
+

γ
n

+

θ
n


)



-


c

2

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α












{



z

12

n



cos





α

+


e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c
n



cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-


a
1


sin





β









1



]








(
113
)













where ζ


n


is substituted by,






ζ


n


=−γ


n


−θ


n








because of the link mechanism to support wheel at this geometric relation.




2.3 Description of the stabilizer linkage point in local coordinate system {1f}.




The stabilizer works as a spring in which force is proportional to the difference of displacement between both arms in a local coordinate system {1f} fixed to the body.













P

stab
.
n


1

f


=






T

4

n


3

n



3

n


1

f




T

8

n


4

n




T

9

n


8

n




T

14

n


9

n




TP

stab
.
n


14

n









=







[



1


0


0


0




0


1


0



b

2

n






0


0


1


0




0


0


0


1



]



[



1


0


0


0




0



cos






γ
n






-
sin







γ
n




0




0



sin






γ
n





cos






γ
n




0




0


0


0


1



]




[



1


0


0


0




0


1


0


0




0


0


1



c

2

n






0


0


0


1



]















[



1


0


0


0




0



cos






θ
n






-
sin







θ
n




0




0



sin






θ
n





cos






θ
n




0




0


0


0


1



]



[



1


0


0


0




0


1


0



e

0

n






0


0


1


0




0


0


0


1



]




[









0




0







0







1



]








=





[



0







e

0

n




cos


(


γ
n

+

θ
n


)



-


c

2

n



sin






γ
n


+

b

2

n










e

0

n




sin


(


γ
n

+

θ
n


)



+


c

2

n



cos






γ
n







0



]








(
114
)













3. Kinetic energy, potential energy and dissipative functions for the <Body>, <Suspension>, <Arm>, <Wheel> and <Stabilizer>.




Kinetic energy and potential energy except by springs are calculated based on the displacement referred to the inertial global coordinate {r}. Potential energy by springs and dissipative functions are calculated based on the movement in each local coordinate.




<Body>










T
b
tr

=


1
2




m
b



(



x
.

b
2

+


y
.

b
2

+


z
.

b
2


)







(
115
)













where
















x
b

=



(


a
0

+

a

1

n



)


cos





β

+


(



b
0


sin





α

+


c
0


cos





α


)


sin





β









y
b

=



b
0


cos





α

-


c
0


sin





α












z
b

=



-

(


a
0

+

a

1

n



)



sin





β

+


(



b
0


sin





α

+


c
0


cos





α


)


cos





β









(
116
)













and





























q

j
,
k


=
β

,
α










x
b




β


=



-

(


a
0

+

a
1


)



sin





β

+


(



b
0


sin





α

+


c
0


cos





α


)


cos





β















x
b




α


=


(



b
0


cos





α

-


c
0


sin





α


)


sin





β














y
b




β


=
0













y
b




α


=



-

b
0



sin





α

-


c
0


cos





α















z
b




β


=



-

(


a
0

+

a
1


)



cos





β

-


(



b
0


sin





α

+


c
0


cos





α


)


sin





β















z
b




α


=


(



b
0


cos





α

-


c
0


sin





α


)


cos





β








(
117
)













and thus













T
b
tr

=






1
2




m
b



(



x
.

b
2

+


y
.

b
2

+


z
.

b
2


)









=






1
2



m
b






j
,
k




(






x
b





q
j








x
b





q
k






q
.

j




q
.

k


+





y
b





q
j








y
b





q
k






q
.

j




q
.

k


+





z
b





q
j








z
b





q
k






q
.

j




q
.

k



)












=






1
2



m
b








β
.

2




{



-

(


a
0

+

a
1


)



sin





β

+


(



b
0


sin





α

+


c
0


cos





α


)


cos





β


}

2


+














α
.

2




{


(



b
0


cos





α

-


c
0


sin





α


)


sin





β

}

2


+
















α
.

2



(



-

b
0



sin





α

-


c
0


cos





α


)


2

+















β
.

2




{



-

(


a
0

+

a
1


)



cos





β

-


(



b
0


sin





α

+


c
0


cos





α


)


sin





β


}

2


+



α
.

2




{


(



b
0


cos





α

-


c
0


sin





α


)


cos





β

}

2


+












2


α
.



β
[



{



-

(


a
0

+

a
1


)



sin





β

+


(



b
0


sin





α

+


c
0


cos





α


)


cos





β


}



(



b
0


cos





α

-


c
0


sin





α


)


sin





β

+

















{



-

(


a
0

+

a
1


)



cos





β

-


(



b
0


sin





α

+


c
0


cos





α


)


sin





β


}



(



b
0


cos





α

-


c
0


sin





α


)


cos





β

]












=






1
2



m
b









α
.

2



(


b
0
2

+

c
0
2


)



+



β
.

2



{



(


a
0

+

a
1


)

2

+


(



b
0


sin





α

+


c
0


cos





α


)

2


}














-
2



α
.




β
.



(


a
0

+

a
1


)




(



b
0


cos





α

-


c
0


sin





α


)










(
118
)










T
b
ro

=






1
2



(



I
bx



ω
bx
2


+


I
by



ω
by
2


+


I
bz



ω
bz
2



)












where












ω
bx

=

α
.














ω
by

=

β
.














ω
bz

=
0


















T
b
ro

=






1
2



(



I
bx




α
.

2


+


I
by




β
.

2



)









U
b

=






m
b



gz
b








=






m
b


g


{



-

(


a
0

+

a
1


)



sin





β

+


(



b
0


sin





α

+


c
0


cos





α


)


cos





β


}









(
119
)









Suspension








T
sn
tr

=


1
2




m
sn



(



x
.

sn
2

+


y
.

sn
2

+


z
.

sn
2


)









where







x
sn

=



{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


sin





β

+


a

1

n



cos





β










y
sn

=



-

z

6

n





sin


(

α
+

γ
n

+

η
n


)



-


c

1

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α







(
120
)








z
sn

=



{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-


a

1

n



sin





β











q

j
,
k


=

z

6

n



,

η
n

,
α
,
β











x
sn





z

6

n




=


cos


(

α
+

γ
n

+

η
n


)



sin





β












x
sn





η
n



=


-

z

6

n





sin


(

α
+

γ
n

+

η
n


)



sin





β












x
sn




α


=


{



-

z

6

n





sin


(

α
+

γ
n

+

η
n


)



-


c

1

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α


}


sin





β












x
sn




β


=



{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-


a

1

n



sin





β













y
sn





z

6

n




=

-

sin


(

α
+

γ
n

+

η
n


)














y
sn





η
n



=


-

z

6

n





cos


(

α
+

γ
n

+

η
n


)














y
sn




α


=



-

z

6

n





cos


(

α
+

γ
n

+

η
n


)



-


c

1

n




cos


(

α
+

γ
n


)



-


b

2

n



sin





α













y
sn




β


=
0





(
121
)











z
sn





z

6

n




=


cos


(

α
+

γ
n

+

η
n


)



cos





β












z
sn





η
n



=


-

z

6

n





sin


(

α
+

γ
n

+

η
n


)



cos





β












z
sn




α


=


{



-

z

6

n





sin


(

α
+

γ
n

+

η
n


)



-


c

1

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α


}


cos





β












z
sn




β


=



-

{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}



sin





β

-


a

1

n



cos





β







(
122
)
















T
sn
tr

=






1
2




m
sn



(



x
.

sn
2

+


y
.

sn
2

+


z
.

sn
2


)









=






1
2



m
sn






j
,
k




(






x
sn





q
j








x
sn





q
k






q
.

j




q
.

k


+





y
sn





q
j








y
sn





q
k






q
.

j




q
.

k


+





z
sn





q
j








z
sn





q
k






q
.

j




q
.

k



)











(
123
)












=






1
2



m
sn








z
.


6

n

2


+



η
.

n
2



z

6

n

2


+



α
.

2

[


z

6

n

2

+

c

1

n

2

+

b

2

n

2

+















2


{



z

6

n




c

1

n



cos






η
n


-


z

6

n




b

2

n




sin


(


γ
n

+

η
n


)



-


c

1

n




b

2

n



sin






γ
n



}


]


+














β
.

2



[



{

(



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


)

}

2

+

a

1

n

2


]


+













2



z
.


6

n




α
.



{



c

1

n



sin






η
n


+


b

2

n




cos


(


γ
n

+

η
n


)




}


-













2



z
.


6

n




β
.







a

1

n




cos


(

α
+

γ
n

+

η
n


)



+













2



η
.

n



α
.



z

6

n




{


z

6

n


+


c

1

n



cos






η
n


-


b

2

n




sin


(


γ
n

+

η
n


)




}


+













2



η
.

n



β
.







z

6

n




a

1

n




sin


(

α
+

γ
n

+

η
n


)



+












2


α
.



β
.







a

1

n




{



z

6

n




sin


(

α
+

γ
n

+

η
2


)



+


c

1

n




sin


(

α
+

γ
n


)



-


b

2

n



cos





α


}










(
124
)










T
sn
ro






0







U
sn

=







m
sn



gz
sn


+


1
2





k
sn



(


z

6

n


-

l
sn


)


2









=







m
sn




gz
sn



[



{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-


a

1

n



sin





β


]



+













1
2





k
sn



(


z

6

n


-

l
sn


)


2









F
sn

=






-

1
2




c
sn




z
.


6

n

2









(
125
)








Arm









T
an
tr

=


1
2




m
an



(



x
.

an
2

+


y
.

an
2

+


z
.

an
2


)








(
126
)






where







x
an

=



{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b
n


sin





α


}


sin





β

+


a

1

n



cos





β










y
an

=



e

1

n




cos


(

α
+

γ
n

+

θ
n


)



-


c

2

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α










z
an

=



{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-


a

1

n



sin





β







(
127
)






and








q

j
,
k


=

θ
n


,
α
,
β











x
an





θ
n



=


e

1

n




cos


(

α
+

γ
n

+

θ
n


)



sin





β












x
an




α


=


{



e

1

n




cos


(

α
+

γ
n

+

θ
n


)



-


c

2

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α


}


sin





β













x
an




β




{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-


a

1

n



sin





β












y
an





θ
n



=


-

e

1

n





sin


(

α
+

γ
n

+

θ
n


)














y
an




α


=



-

e

1

n





sin


(

α
+

γ
n

+

θ
n


)



-


c

2

n




cos


(

α
+

γ
n


)



-


b

2

n



sin





α













y
an




β


=
0











z
an





θ
n



=


e

1

n




cos


(

α
+

γ
n

+

θ
n


)



cos





β












z
an




α


=


{



e

1

n




cos


(

α
+

γ
n

+

θ
n


)



-


c

2

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α


}


cos





β












z
an




β


=



-

{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}



sin





β

-


a

1

n



cos





β








thus




(
128
)










T
an
tr

=






1
2




m
an



(



x
.

an
2

+


y
.

an
2

+


z
.

an
2


)









=






1
2



m
an






j
,
k




(






x
an





q
j








x
an





q
k






q
.

j




q
.

k


+





y
an





q
j








y
an





q
k






q
.

j




q
.

k


+





z
an





q
j








z
an





q
k






q
.

j




q
.

k



)










(
129
)









=






1
2



m
an








θ
.

n
2



e

1

n

2


+



α
.

2

[


e

1

n

2

+

c

2

n

2

+

b

2

n

2

-

2


{



e

1

n




c

2

n



sin






θ
n


+


e

1

n




b

2

n




cos


(


γ
n

+

θ
n


)



+























c

2

n




b

2

n



sin






γ
n


}


]

+



β
.

2



[



{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}

2

+

a

1

n

2


]


+












2


θ
.



α
.



e

1

n




{


e

1

n


-


c

2

n



sin






θ
n


+


b

2

n




cos


(


γ
n

+

θ
n


)




}


-

2



θ
.

n



β
.



e

1

n




a

1

n




cos


(

α
+

γ
n

+

θ
n


)



-













2

αβ






a

1

n




{



e

1

n




cos


(

α
+

γ
n

+

θ
n


)



-


c

1

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α


}











(
130
)










T
an
ro

=






1
2



I
ax



ω
ax
2








=






1
2





I
ax



(


α
.

+


θ
.

n


)


2









U
an

=






m
an



gz
an








=






m
an



g


[



{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-


a

1

n



sin





β


]










(
131
)








Wheel









T
wn
tr

=


1
2




m
wn



(



x
.

wn
2

+


y
.

wn
2

+


z
.

wn
2


)








(
132
)













where
















x
wn

=



{



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


sin





β

+


a

1

n



cos





β









y
wn

=



e

3

n




cos


(

α
+

γ
n

+

θ
n


)



-


c

2

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α












z
wn

=



{



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-


a

1

n



sin





β









(
133
)













Substituting m


an


with m


wn


and e


1n


with e


3n


in the equation for arm, yields an equation for the wheel as:













T
wn
tr

=


1
2



m
wn








θ
.

n
2



e

3

n

2


+



α
.

2



[


e

3

n

2

+

c

2

n

2

+

b

2

n

2

-

2


{



e

2

n




c

2

n



sin






θ
n


+


e

3

n




b

2

n




cos


(


γ
n

+

θ
n


)



+


c

2

n




b

2

n



sin






γ
n



}



]


+



β
.

2



[



{



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}

2

+

a

1

n

2


]


+

2


θ
.



α
.



e

3

n




{


e

3

n


-


c

2

n



sin






θ
n


+


b

2

n




cos


(


γ
n

+

θ
n


)




}


-

2



θ
.

n



β
.



e

3

n




a

1

n




cos


(

α
+

γ
n

+

θ
n


)



-

2


α
.



β
.



a

1

n




{



e

3

n




cos


(

α
+

γ
n

+

θ
n


)



-


c

1

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α


}







(
134
)










T
wn
ro

=




0







U
wn

=







m
wn



gz
wn


+


1
2





k
wn



(


z

12

n


-

l
wn


)


2









=






m
wn



g
[



{



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-

















a

1

n



sin





β

]


+


1
2





k
wn



(


z

12

n


-

l
wn


)


2









F
wn

=






-

1
2




c
wn




z
.


12

n

2









(
135
)








Stabilizer









T
zn
tr


0





(
136
)







T
zn
ro


0




(
137
)











U
zn








1
2





k
zn



(


z
zi

-

z
zii


)


2








=






1
2



k
zn




{



e

0

i




sin


(


γ
i

+

θ
i


)



-


e

0

ii




sin


(


γ
ii

+

θ
ii


)




}

2








=






1
2



k
zn



e

0

i

2




{


sin


(


γ
i

+

θ
i


)


+

sin


(


γ
i

+

θ
i


)



}

2











where







e

0

ii


=

-

e

0

i








(
138
)







F
zn


0




(
139
)













Therefore the total kinetic energy is:










T
tot

=






T
b
tr

+



n

i
,
ii




&LeftBracketingBar;


T
sn
tr

+

T
an
tr

+

T
wn
tr

+

T
b
ro

+

T
an
ro









(
140
)










T
tot

=






T
b
tr

+



n

i
,
ii




&LeftBracketingBar;


T
sn
tr

+

T
an
tr

+

T
wn
tr

+

T
b
ro

+

T
an
ro


&RightBracketingBar;












=






1
2



m
b









α
.

2



(


b
0
2

+

c
0
2


)



+














β
.

2



{



(


a
0

+

a
1


)

2

+


(



b
0


sin





α

+


c
0


cos





α


)

2


}


-













2


α
.




β
.



(


a
0

+

a
1


)




(



b
0


cos





α

-


c
0


sin





α


)





+


n

i
,
ii



&RightBracketingBar;



1
2



m
sn




z
.


6

n

2


+















η
.

n
2



z

6

n

2


+



α
.

2

[

z

6

n

2


+

c

1

n

2

+

b

2

n

2

+

2


{



z

6

n




c

1

n



cos






η
n


-























z

6

n




b

2

n




sin


(


γ
n

+

η
n


)



-


c

1

n




b

2

n



sin






γ
n



}


]

+


c

1

n




b

2

n



sin






γ
n



}

]

+













β
.

2

[

{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+





















c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


2

+

a

1

n

2


]

+












2



z
.


6

n




α
.



{



c

1

n



sin






η
n


+


b

2

n




cos


(


γ
n

+

η
n


)




}


-













2







z
.


6

n




β
.







a

1

n




cos


(

α
+

γ
n

+

η
n


)



+













2



η
.

n



α
.







z

6

n




{


z

6

n


+


c

1

n



cos






η
n


-


b

2

n




sin


(


γ
n

+

η
n


)




}


+













2



η
.

n



β
.







z

6

n




a

1

n




sin


(

α
+

γ
n

+

η
n


)



+












2


α
.



β
.



a

1

n




{



z

6

n




sin


(

α
+

γ
n

+

η
n


)



+


c

1

n




sin


(

α
+

γ
n


)



-


















b

2

n



cos





α

}







+

1
2




m
an







θ
.

n
2



e

1

n

2


+



α
.

2

[


e

1

n

2

+

c

2

n

2

+

b

2

n

2

-













2


{



e

1

n




c

2

n



sin






θ
n


-


e

1

n




b

2

n




cos


(


γ
n

+

θ
n


)



+


















c

2

n




b

2

n



sin






γ
n


}


]

+



β
.

2

[

{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+





















c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


2



a

1

n

2


]

+












2


θ
.







α
.



e

1

n




{


e

1

n


-


c

2

n



sin






θ
n


+


b

2

n




cos


(


γ
n

+

θ
n


)




}


-













2



θ
.

n



β
.



e

1

n




a

1

n




cos


(

α
+

γ
n

+

θ
n


)



-












2






α
.



β
.



a

1

n




{



e

1

n




cos


(

α
+

γ
n

+

θ
n


)



-
















c

1

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α





+

1
2




m
wn







θ
.

n
2



e

3

n

2


+















α
.

2

[

e

3

n

2


+

c

2

n

2

+

b

2

n

2

-

2


{



e

3

n




c

2

n



sin






θ
n


-




















e

3

n




b

2

n




cos


(


γ
n

+

θ
n


)



+


c

2

n




b

2

n



sin






γ
n



}


]

+













β
.

2

[

{



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+




















b

2

n



sin





α

}


2

+

a

1

n

2


]

+












2


θ
.



α
.



e

3

n




{


e

3

n


-


c

2

n



sin






θ
n


+


b

2

n




cos


(


γ
n

+

θ
n


)




}


-













2



θ
.

n



β
.



e

3

n




a

1

n




cos


(

α
+

γ
n

+

θ
n


)



-












2


α
.



β
.



a

1

n




{



e

3

n




cos


(

α
+

γ
n

+

θ
n


)



-


















c

1

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α


}





+















1
2



(



I
bx




α
.

2


+


I
by




β
.

2



)


+


1
2





I
anx



(


α
.

+


θ
.

n


)


2



&RightBracketingBar;








(
141
)









     


=






1
2

[




α
.

2



m
bbI


+



β
.

2



{


m
baI

+



m
b



(



b
0


sin





α

+


c
0


cos





α


)


2


}


-
















2


α
.



β
.




m
ba



(



b
0


cos





α

-


c
0


sin





α


)



]


+












1
2






n


i
,
ii




&LeftBracketingBar;



m
sn



(



z
.


6

n

2

+



η
.

n
2



z

6

n

2



)


+



θ
.

n
2



m
aw2In


+



















α
.

2






m
sawIn


+


m
sn




z

6

n


[


z

6

n


+

2


m
sn



{



c

1

n



cos






η
n


-





















b

2

n




sin


(


γ
n

+

η
n


)



}


]

-

2


m
aw1n



{



c

2

n



sin






θ
n


-

















b

2

n




cos


(


γ
n

+

θ
n


)



}








+


β
.

2




m
saw2n


+














m
sn



{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+



















c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


2

+













m
an




{



e
1



sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}

2


+













m
wn




{



e
3



sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}

2




+














2



z
.


6

n




α
.



m
sn



{



c

1

n



sin






η
n


+


b

2

n




cos


(


γ
n

+

η
n


)




}


-













2



z
.


6

n




β
.







ma

1

n




cos


(

α
+

γ
n

+

η
n


)



+













2



η
.

n



α
.







m
sn



z

6

n




{


z

6

n


+


c

1

n



cos






η
n


-


b

2

n




sin


(


γ
n

+

η
n


)




}


+













2



η
.

n



β
.







m
sn



z

6

n




a

1

n




sin


(

α
+

γ
n

+

η
n


)



+













2


θ
.








α
.



[


m
aw2In

-


m
aw1n



{



c

2

n



sin






θ
n


-


b

2

n




cos


(


γ
n

+

θ
n


)




}



]



-













2


θ
.







β
.



m
aw1n



a

1

n




cos


(

α
+

γ
n

+

θ
n


)



+












2


α
.



β
.



m
aw1n



{



m
sawcn



sin


(

α
+

γ
n


)



-


m
sawbn


cos





α

+


















m
sn



z

6

n




sin


(

α
+

γ
n

+

η
n


)



-


m
aw1n



cos


(

α
+

γ
n

+

θ
n


)




}


&RightBracketingBar;







(
142
)








where









m
ba

=






m
b



(


a
0

+

a
1


)









m
bbI

=







m
b



(


b
0
2

+

c
0
2


)


+

I
bx









m
baI

=








m
b



(


a
0

+

a
1


)


2

+

I
by









m
sawan

=






(


m
sn

+

m
an

+

m
wn


)



a

1

n










m
sawbn

=






(


m
sn

+

m
an

+

m
wn


)



b

2

n










m
sawcn

=







m
sn



c

1

n



+


(


m
an

+

m
wn


)



c

2

n











m
saw2n

=






(


m
sn

+

m
an

+

m
wn


)



a

1

n

2









m
sawIn

=







m
an



e

1

n

2


+


m
wn



e

3

n

2


+


m
sn



(


c

1

n

2

+

b

2

n

2

-

2


c

1

n




b

2

n



sin






γ
n



)


+














(


m
an

+

m
wn


)



(


c

2

n

2

+

b

2

n

2

-

2


c

2

n




b

2

n



sin






γ
n



)


+

I
axn









m
aw2In

=







m
an



e

1

n

2


+


m
wn



e

3

n

2


+

I
axn









m
aw1n

=







m
an



e

1

n



+


m
wn



e

3

n











m
aw2n

=







m
an



e

1

n

2


+


m
wn



e

3

n

2













(
143
)













Hereafter variables and constants which have index ′n′ implies implicit or explicit that they require summation with n=i, ii.




Total potential energy is:










U
tot

=



U
b

+


n

i
,
ii



|


U
sn

+

U
an

+

U
wn

+

U
zn







(
144
)









=







m
b


g


{



-

(


a
0

+

a
1


)



sin





β

+


(



b
0


sin





α

+


c
0


cos





α


)


cos





β


}


+














n

i
,
ii




&LeftBracketingBar;


m
sn


g
[

{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+






















b

2

n



sin





α

}



cos





β

-


a

1

n



sin





β


]

+


1
2





k
sn



(


z

6

n


-

l
sn


)


2


+












m
an


g
[

{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+




















b

2

n



sin





α

}



cos





β

-


a

1

n



sin





β


]

+












m
wn


g
[

{



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+




















b

2

n



sin





α

}



cos





β

-


a

1

n



sin





β


]

+


1
2





k
wn



(


z

12

n


-

l
wn


)


2


+












1
2



k
zn

[



e

0

n




{


sin


(


γ
i

+

θ
i


)


+

sin


(


γ
ii

+

θ
ii


)



}


+

















c

2

n




(


cos






γ
i


-

cos






γ
ii



)


]


2

&RightBracketingBar;







(
145
)










=






g


{



-

m
ba



sin





β

+



m
b



(



b
0


sin





α

+


c
0


cos





α


)



cos





β


}


+














1
2



k
zi



e

0

i

2




{


sin


(


γ
i

+

θ
i


)


+

sin


(


γ
ii

+

θ
ii


)



}

2


+














n

i
,
ii






g
[

{



m
sn



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+
















m
aw1n



sin


(

α
+

γ
n

+

θ
n


)



+


m
sawcn



cos


(

α
+

γ
n


)



+


















m
sawbn


sin





α

}



cos





β

-


m
sawan


sin





β


]

+













1
2





k
sn



(


z

6

n


-

l
sn


)


2


+


1
2





k
wn



(


z

12

n


-

l
wn


)


2













where




(
146
)



















m
ba

=


m
b



(


a
0

+

a
1


)









m
sawan

=


(


m
sn

+

m
an

+

m
wn


)



a

1

n













m
sawbn

=


(


m
sn

+

m
an

+

m
wn


)



b

2

n













m
sawcn

=



m
sn



c

1

n



+


(


m
an

+

m
wn


)



c

2

n














γ
ii

=

-

γ
i









(
147
)













4. Lagrangye's Equation




The Lagrangian is written as:












L
=






T
tot

-

U
tot








=






1
2

[




α
.

2



m
bbI


+



β
.

2



{


m
baI

+



m
b



(



b
0


sin





α

+


c
0


cos





α


)


2


}


-















2


α
.



β
.




m
ba



(



b
0


cos





α

-


c
0


sin





α


)



]


+


1
2





n

i
,
ii




&LeftBracketingBar;



m
sn



(



z
.


6

n

2

+



η
.

n
2



z

6

n

2



)


+





















θ
.

n
2



m
aw2In


+


α
.

2







m
sawIn


+


m
sn




z

6

n


[


z

6

n


+

2


{



c

1

n



cos






η
n


-





















b

2

n




sin


(


γ
n

+

η
n


)



}


]

-

2


m
aw1n



{



c

2

n



sin






θ
n


-


b

2

n



cos






(

γ
+

θ
n


)



}




+

















β
.

2






m
saw2n


+


m
sn



{



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+


c

1

n




cos


(

α
+

γ
n


)



+


















b

2

n



sin





α

}


2

+


m

an








{



e

1

n



sin






(

α
+

γ
n

+

θ
n


)


+


c

2

n



cos






(

α
+

γ
n


)


+


















b

2

n



sin





α

}


2

+


m

2

n




{



e
3


sin






(

α
+

γ
n

+

θ
n


)


+


c

2

n



cos






(

α
+

γ
n


)


+


















b

2

n



sin





α

}


2







+
2




z
.


6

n




α
.



m
sn



{



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+

η
n


)



}


-














2



z
.


6

n




β
.



m
sn



a

1

n



cos






(

α
+

γ
n

+

η
n


)


+

2



η
.

n



α
.



m
sn



z

6

n




{


z

6

n


+



















c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+

2



η
.

n



β
.



m
sn



z

6

n




a

1

n



sin






(

α
+

















γ
n

+

η
n


)


+

2


θ
.




α
.



[


m
aw2In



{



c

2

n



sin






θ
n


-


b

2

n




cos


(


γ
n

+

θ
n


)




}


]



-












2


θ
.



β
.



m
aw1n



a

1

n



cos






(

α
+

γ
n

+

θ
n


)


+

2


α
.



β
.



a

1

n




{


m
sawcn



sin
(

α
+



















γ
n

)


-


m
sawbn


cos





α

+


m
sn



z

6

n



sin






(

α
+

γ
n

+

η
n


)


-















m
aw1n


cos






(

α
+

γ
n

+

θ
n


)


}


&RightBracketingBar;

-

g


{



-

m
ba



sin





β

+


m
b

(



b
0


sin





α

+





















c
0


cos





α

)



cos





β

}

+


1
2



k
zi



e

0

i

2




{


sin






(


γ
i

+

θ
i


)


+

sin






(


γ
ii

+

θ
ii


)



}

2


-













n

i
,
ii






g
[

{



m
sn



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+


m
aw1n


sin






(

α
+

γ
n

+

θ
n


)


+





















m
sawcn


cos






(

α
+

γ
n


)


+


m
sawbn


sin





α


}



cos





β

-


m
sawan


sin





β


]

+













1
2





k
sn



(


z

6

n


-

l
sn


)


2


+


1
2





k
wn



(


z

12

n


-

l
wn


)


2











(
148
)












L



β


=









g


{



m
ba


cos





β

+


m
b



{



b
0


sin





α

+


c
0


cos





α







)


sin





β

}

+














n

i
,
ii






g
[

{



m
sn



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+
















m
aw1n


sin






(

α
+

γ
n

+

θ
n


)


+


m
sawcn


cos






(

α
+

γ
n


)


+


















m
sawbn


sin





α

}



sin





β

+


m
sawan


cos





β


]









(
149
)

























L



α


=






{




β
.

2




m
b



(



b
0


cos





α

-


c
0


sin





α


)



+


α
.



β
.



m
ba



}



(



b
0


sin





α

+



















c
0


cos





α

)


+



n

i
,
ii




&LeftBracketingBar;





+


β
.

2









m
sn



{



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+

















c


1

n








cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}




{



-

z

6

n




sin






(

α
+

γ
n

+

η
n


)


-

















c

1

n



sin






(

α
+

γ
n


)


+


b

2

n



cos





α


}


+


m
an



{



e

1

n



sin






(

α
+

γ
n

+

θ
n


)


+


















c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}




{



e
1


cos






(

α
+

γ
n

+

θ
n


)


-

















c

2

n



sin






(

α
+

γ
n


)


+


b

2

n



cos





α


}


+


m

2

n




{



e

3

n



sin






(

α
+

γ
n

+

θ
n


)


+


















c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}




{



e
3


cos






(

α
+

γ
n

+

θ
n


)


-

















c

2

n



sin






(

α
+

γ
n


)


+


b

2

n



cos





α


}







+


z
.


6

n





β
.



m
sn



a

1

n



sin






(

α
+

γ
n

+

















η
n

)


+



η
.

n



β
.



m
sn



z

6

n




a

1

n



cos






(

α
+

γ
n

+

η
n


)


+













θ
.



β
.



m
aw1n



a

1

n



sin






(

α
+

γ
n

+

θ
n


)


+


α
.



β
.



a

1

n




{


m
sawcn


cos






(

α
+



















γ
n

)


+


m
sawbn


sin





α

+


m
sn



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+















m
aw1n


sin






(

α
+

γ
n

+

θ
n


)


}


&RightBracketingBar;

-


gm
b

(



b
0


cos





α

-

















c
0


sin





α

)



cos





β

+



n

i
,
ii




g
(



m
sn



z

6

n



sin






(

α
+

γ
n

+

η
n


)


-
















m
aw1n


cos






(

α
+

γ
n

+

θ
n


)


+


m
sawcn


sin






(

α
+

γ
n


)


-















m
sawbn


cos





α

}



cos





β







(
150
)

























L




η
n



=








α
.

2



m
sn



z

6

n




{



-

c

1

n




sin






η
n


-


b

2

n



cos






(


γ
n

+

η
n


)



}


+














β
.

2



m
sn



{



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+


c

1

n



cos






(

α
+

γ
n


)


+


















b

2

n



sin





α

}




{


-

z

6

n




sin






(

α
+

γ
n

+

η
n


)


}


+



z
.


6

n




α
.



m
sn



{



c

1

n



cos





η

-

















b

2

n



sin






(


γ
n

+

η
n


)


}


+



z
.


6

n




β
.



m
sn



a

1

n



sin






(

α
+

γ
n

+

η
n


)


-














η
.

n



α
.



m
sn



z

6

n




{



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+

η
n


)



}


+















η
.

n



β
.



m
sn



z

6

n




a

1

n



cos






(

α
+

γ
n

+

η
n


)


+














α
.



β
.



a

1

n




m
sn



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+
g













gm
sn



z

6

n



sin






(

α
+

γ
n

+

η
n


)


cos





β








(
151
)












L




θ
n



=







-

k
zi




e

0

i

2



{


sin






(


γ
i

+

θ
i


)


+

sin






(


γ
ii

+

θ
ii


)



}


cos






(


γ
n

+

θ
n


)


-















α
.

2



m
aw1n



{



c

2

n



cos






θ
n


+


b

2

n



sin






(


γ
n

+

θ
n


)



}


+














β
.

2






m
an



{



e

1

n



sin






(

α
+

γ
n

+

θ
n


)


+


c

2

n



cos






(

α
+

γ
n


)


+





















b

2

n



sin





α

}



1

n



cos






(

α
+

γ
n

+

θ
n


)


+


m
wn



{


e

3

n



sin






(

α
+





















γ
n

+

θ
n


)


+


c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}



e

3

n



cos






(

α
+

















γ
n

+

θ
n


)




-


θ
.



α
.



m
aw1n



{



c

2

n



cos






θ
n


+


b

2

n



sin






(


γ
n

+

θ
n


)



}


+













θ
.



β
.



m
aw1n



a

1

n



sin






(

α
+

γ
n

+

θ
n


)


+


α
.



β
.



a

1

n




m
aw1n


sin






(

α
+


















γ
n

+

θ
n


)


-


gm
aw1n


cos






(

α
+

γ
n

+


θ
n


cos





β











(
152
)

























L




z

6

n




=







m
sn




η
.

n
2



z

6

n



+



α
.

2




m
sn

[


z

6

n


+

{



c


1

n








cos






η
n


-





















b

2

n



sin






(


γ
n

+

η
n


)


}


]

+



β
.

2



m
sn



{



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+



















c

1

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}



cos






(

α
+

γ
n

+

η
n


)


+














η
.

n



α
.



m
sn



{


2


z

6

n



+


c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+















η
.

n



β
.



m
sn



a

1

n



sin






(

α
+

γ
n

+

η
n


)


+


α
.



β
.



a

1

n




m
sn


sin






(

α
+


















γ
n

+

η
n


)


-


gm
sn


cos






(

α
+

γ
n

+

η
n


)


cos





β

-












k
sn



(


z

6

n


-

l
sn


)









(
153
)









L




z

12

n




=

-


k
wn



(


z

12

n


-

l
wn


)







(
154
)












L




β
.



=






β
.






m
saw2n

+

m
baI

+



m
b



(



b

0







sin





α

+


c
0


cos





α


)


2

+















m
sn



{



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+


c

1

n



cos






(

α
+

γ
n


)


+


















b

2

n



sin





α

}


2

+


m
an



{



e

1

n



sin






(

α
+

γ
n

+

θ
n


)


+



















c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}


2

+


m
wn



{


e

3

n



sin






(

α
+
















(









γ
n

+

θ
n


)


+


c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α



}

)

2



-













α
.




m
ba



(



b
0


cos





α

-


c
0


sin





α


)



-



z
.


6

n




m
sn



a

1

n



cos






(

α
+


















γ
n

+

η
n


)


+



η
.

n



m
sn



z

6

n




a

1

n



sin






(

α
+

γ
n

+

η
n


)


-













θ
.



m
aw1n



a

1

n



cos






(

α
+

γ
n

+

θ
n


)


+













α
.



a

1

n




{



m
sawcn


sin






(

α
+

γ
n


)


-


m
sawbn


cos





α

+

















m
sn



z

6

n



sin






(

α
+

γ
n

+

η
n


)


-


m
aw1n


cos






(

α
+

γ
n

+

θ
n


)



}








(
155
)






























t




(



L




β
.



)


=





β
¨






m
saw2n


+

m
baI

+



m
b



(



b
0


sin





α

+


c
0


cos





α


)


2

+












m
sn



{



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+


c

1

n



cos






(

α
+

γ
n


)


+


















b

2

n



sin





α

}


2

+


m
an



{



e

1

n



sin






(

α
+

γ
n

+

θ
n


)


+



















c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}


2

+


m

2

n




{


e

3

n



sin






(

α
+















(









γ
n

+

θ
n


)


+


c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α



}

)

2




+















2


β
.







α
.




m
b



(



b
0


sin





α

+


c
0


cos





α


)




(



b
0


cos





α

-


c
0


sin





α


)


+













m
.

sn



{



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+


c

1

n



cos






(

α
+

γ
n


)


+

















b

2

n



sin





α

}




{




z
.


6

n



cos






(

α
+

γ
n

+

η
n


)


-















(


α
.

+


η
.

n


)



z

6

n



sin






(

α
+

γ
n

+

η
n


)


-


α
.

[



c

1

n



sin






(

α
+

γ
n


)


-


















b

2

n



cos





α

]


}

+


m
an



{



e

1

n



sin






(

α
+

γ
n

+

θ
n


)


+


















c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}




{


(


α
.

+


θ
.

n


)



e

1

n



cos






(

α
+




















γ
n

+

θ
n


)


-


α
.



[



c

2

n



sin






(

α
+

γ
n


)


-


b

2

n



cos





α


]



}

+












m
wn



{



e

3

n



sin






(

α
+

γ
n

+

θ
n


)


+


c

2

n



cos






(

α
+

γ
n


)


+

















b

2

n



sin





α

}




{



(


α
.

+


θ
.

n


)



e

3

n



sin






(

α
+

γ
n

+

θ
n


)


-
















α
.



[



c

2

n



sin






(

α
+

γ
n


)


-


b

2

n



cos





α


]


}





-














α
¨




m
ba



(



b
0


cos





α

-


c
0


sin





α


)



+



α
.

2




m
ba

(



b
0


sin





α

+


















c
0


cos





α

)


-



z
¨


6

n




m
sn



a

1

n



cos






(

α
+

γ
n

+

η
n


)


+















z
.


6

n




(


α
.

+


η
.

n


)




m
sn



a

1

n



sin






(

α
+

γ
n

+

η
n


)


+















η
¨

n



m
sn



z

6

n




a

1

n



sin






(

α
+

γ
n

+

η
n


)


+















η
.

n



m
sn




z
.


6

n




a

1

n



sin






(

α
+

γ
n

+

η
n


)


+



η
.

n

(


α
.

+


















η
.

n

)




m
sn



z

6

n




a

1

n



cos






(

α
+

γ
n

+

η
n


)


-














θ
¨

n



m
aw1n



a

1

n



cos






(

α
+

γ
n

+

θ
n


)


+



θ
.

n

(


α
.

+


















θ
.

n

)




m
aw1n



a

1

n



sin






(

α
+

γ
n

+

θ
n


)


+












α
¨



a

1

n




{



m
sawcn


sin






(

α
+

γ
n


)


-


m
sawbn


cos





α

+
















m
sn



z

6

n



sin






(

α
+

γ
n

+

η
n


)


-


m
aw1n


cos






(

α
+

γ
n

+


















θ
n

)


}

+


α
.



a

1

n




{



α
.



m
sawcn


cos






(

α
+

γ
n


)


+


α
.



m
sawbn


sin





α

+
















(


α
.

+


η
.

n


)



m
sn



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+














m
sn




z
.


6

n



sin






(

α
+

γ
n

+

η
n


)


+


(


α
.

+


θ
.

n


)



m
aw1n


sin






(

α
+


















γ
n

+

θ
n


)


}







(
156
)




























L




α
.



=







α
.



m
bbI


-


β
.




m
ba



(



b
0


cos





α

-


c
0


sin





α


)



+

&AutoLeftMatch;

α


&AutoLeftMatch;
.









m
sawIn


+













m
sn




z

6

n




[


z

6

n


+

2


{



c

1

n



cos





η

-


b

2

n



sin






(


γ
n

+

η
n


)



}



]



-












2


m
aw1n



{



c

2

n



sin






θ
n


-


b

2

n



cos






(


γ
n

+

θ
n


)



}




+
















z
.


6

n




m
sn



{



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+

η
n


)



}


+















η
.

n



m
sn



z

6

n




{


z

6

n


+


c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+














θ
.



[


m
aw2In

-


m
aw1n



{



c

2

n



sin






θ
n


-


b

2

n



cos






(


γ
n

+

θ
n


)



}



]


+













β
.



a

1

n




{



m
sawcn


sin






(

α
+

γ
n


)


-


m
sawbn


cos





α

+

















m
sn



z

6

n



sin






(

α
+

γ
n

+

η
n


)


-


m
aw1n


cos






(

α
+

γ
n

+

θ
n


)



}








(
157
)














t




(



L




α
.



)


=







-

β
¨





m
ba



(



b
0


cos





α

-


c
0


sin





α


)



+


β
.



α
.




m
ba

(



b
0


sin





α

+





















c
0


cos





α

)


+

α
¨






m
bbI


+

m
sawIn

+


m
sn




z

6

n


[


z

6

n


+
















2






{



c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


]


-











2


m
aw1n



{



c

2

n



sin






θ
n


-


b

2

n



cos






(


γ
n

+

θ




)






}




+















α
.






m
sn





z
.


6

n


[


z

6

n


+

2


{



c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+























η
n

)


}

]

+


m
sn




z

6

n


[



z
.


6

n


-

2



η
.

n



{



c

1

n



sin






η
n


+





















b

2

n



cos






(


γ
n

+

η
n


)


}


]

-

2



θ
.

n



m
aw1n



{



c

2

n



cos






θ
n


+

















b

2

n



sin






(


γ
n

+

θ
n


)


}







+


z
¨


6

n





m
sn



{



c

1

n



sin






η
n


+


















b

2

n



cos






(


γ
n

+

η
n


)


}


+



z
.


6

n





η
.

n



m
sn



{



c

1

n



cos






η
n


-

















b

2

n



sin






(


γ
n

+

η
n


)


}


+



η
¨

n



m
sn



z

6

n




{


z

6

n


+


















c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+














η
.

n



m
sn




z
.


6

n




{


z

6

n


+


c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+














η
.

n



m
sn



z

6

n




{



z
.


6

n


-



η
.

n

[



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+























η
n

)


]

}

+


θ
¨

[


m
aw2In

-


m
aw1n



{



c

2

n



sin






θ
n


-




















b

2

n



cos






(


γ
n

+

θ
n


)


}


]

-



θ
.

n
2



m
aw1n



{



c

2

n



cos






θ
n


+


















b

2

n



sin






(


γ
n

+

θ
n


)


}


]

+


β
¨



a

1

n




{



m
sawcn


sin






(

α
+

γ
n


)


-
















m
sawbn


cos





α

+


m
sn



z

6

n



sin






(

α
+

γ
n

+

η
n


)


-















m
aw1n


cos






(

α
+

γ
n

+

θ
n


)


}


+


β
.



a

1

n




{


α
.

[


m
sawcn


cos






(

α
+





















γ
n

)


+


m
sawbn


sin





α


]

+


m
sn




z
.


6

n



sin






(

α
+

γ
n

+

η
n


)


+













(


α
.

+


η
.

n


)



m
sn



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+

(


α
.

+


















θ
.

n

)




m
aw1n


sin






(

α
+

γ
n

+

θ
n


)


}







(
158
)












L





η
.

n



=







m
sn




η
.

n



z

6

n

2


+


α
.



m
sn



z

6

n




{


z

6

n


+


c

1

n



cos






η
n


-


















b

2

n



sin






(


γ
n

+

η
n


)


}


+


β
.



m
sn



z

6

n




a

1

n



sin






(

α
+

γ
n

+

η
n


)









(
159
)



























t




(



L





η
.

n



)


=







m
sn




η
¨

n



z

6

n

2


+

2


m
sn




η
.

n




z
.


6

n




z

6

n



+














α
¨



m
sn



z

6

n




{


z

6

n


+


c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+














α
.



m
sn




z
.


6

n




{


z

6

n


+


c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+














α
.



m
sn



z

6

n




{



z
.


6

n


-



η
.

n



[



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+

η
n


)



]



}


+














β
¨



m
sn



z

6

n




a

1

n



sin






(

α
+

γ
n

+

η
n


)


+














β
.



m
sn




z
.


6

n




a

1

n



sin






(

α
+

γ
n

+

η
n


)


+














β
.



(


α
.

+


η
.

n


)




m
sn



z

6

n




a

1

n



cos






(

α
+

γ
n

+

η
n


)









(
160
)












L





θ
.

n



=








θ
.

n



m
aw2In


+


α
.

[


m
aw2In

-


m
aw1n



{



c

2

n



sin






θ
n


-





















b

2

n



cos






(


γ
n

+

θ
n


)


}


]

-


β
.



m
aw1n



a

1

n



cos






(

α
+

γ
n

+

θ
n


)









(
161
)














t




(



L





θ
.

n



)


=








θ
¨

n



m
aw2In


+


α
¨

[


m
aw2In

-


m
aw1n



{



c

2

n



sin






θ
n


-





















b

2

n



cos






(


γ
n

+

θ
n


)


}


]

-


α
.




θ
.

n



m
aw1n



{



c

2

n



cos






θ
n


+

















b

2

n



sin






(


γ
n

+

θ
n


)


}


-


β
¨



m
aw1n



a

1

n



cos






(

α
+

γ
n

+
















θ
n

)


+



β
.



(


α
.

+


θ
.

n


)




m
aw1n



a

1

n



sin






(

α
+

γ
n

+

θ
n


)









(
162
)












L





z
.


6

n




=







m
sn




z
.


6

n



+


α
.



m
sn



{



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+

η
n


)



}


-













β
.



m
sn



a

1

n



cos






(

α
+

γ
n

+

η
n


)









(
163
)














t




(



L





z
.


6

n




)


=







m
sn




z
¨


6

n



+


α
¨



m
sn



{



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+

η
n


)



}


+














α
.




η
.

n



m
sn



{



c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


-














β
¨



m
sn



a

1

n



cos






(

α
+

γ
n

+

η
n


)


+














β
.



(


α
.

+


η
.

n


)




m
sn



a

1

n



sin






(

α
+

γ
n

+

η
n


)









(
164
)









L





z
.


12

n




=
0




(
165
)











t




(



L





z
.


12

n




)


=
0




(
166
)













The dissipative function is:










F
tot

=


-

1
2




(



c
sn




z
.


6

n

2


+


c
wn




z
.


12

n

2



)






(
167
)













The constraints are based on geometrical constraints, and the touch point of the road and the wheel, The geometrical constraint is expressed as






e


2n


cos θ


n


=−(z


6n


−d


1n


)sin η


n


e


2n


sin θ


n


−(z


6n


−d


1n


)cos η


n


=c


1n


−c


2n


  (168)






The touch point of the road and the wheel is defined as













z
tn

=





z

P

touchpoint





n










=





{



z

12

n



cos





α

+


e

3

n



sin






(

α
+

γ
n

+

θ
n


)


+


c

2

n



cos






(

α
+

γ
n


)


+

















b

2

n



sin





α

}



cos





β

-


a

1

n



sin





β







=






R
n



(
t
)









(
169
)













where R


n


(t) is road input at each wheel.




Differentials are:














θ
.

n



e

2

n



sin






θ
n


-



z
.


6

n



sin






η
n


-




η
.

n



(


z

6

n


-

d

1

n



)



cos






η
n



=
0











θ
.

n



e

2

n



cos






θ
n


-



z
.


6

n



cos






η
n


+




η
.

n



(


z

6

n


-

d

1

n



)



sin






η
n



=
0










{




z
.


12

n



cos





α

-


α
.



z

12

n



sin





α

+


(


α
.

+


θ
.

n


)



e

3

n




cos


(

α
+

γ
n

+

θ
n


)



-


α
.



c

2

n




sin


(

α
+

γ
n


)



+


α
.



b

2

n



cos





α


}


cos





β

-


β
.



[



{



z

12

n



cos





α

+


e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


sin





β

+


a

1

n



cos





β


]


-



R
.

n



(
t
)



=
0





(
170
)













Since the differentials of these constraints are written as













j




a
lnj


d







q
.

j



+


a
lnt


dt


=

0






(


l
=
1

,
2
,


3





n

=
i

,
ii

)






(
171
)













then the values a


lnj


are obtained as follows.












a

1

n1


=
0

,


a

1

n2


=
0

,


a

1

n3


=


-

(


z

6

n


-

d

1

n



)



cos






η
n



,


a

1

n4


=


e

2

n



sin






θ
n



,


a

1

n5


=


-
sin







η
n



,


a

1

n6


=
0










a

2

n1


=
0

,


a

2

n2


=
0

,


a

2

n3


=


(


z

6

n


-

d

1

n



)


sin






η
n



,


a

2

n4


=


e

2

n



cos






θ
n



,


a

2

n5


=


-
cos







η
n



,


a

2

n6


=
0










a

3

n1


=



-

{



z

12

n



cos





α

+


e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}



sin





β

+


a

1

n



cos





β



,






a

3

n2


=


{



-

z

12

n




sin





α

+


e

3

n




cos


(

α
+

γ
n

+

θ
n


)



-


c

2

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α


}


cos





β


,






a

3

n3


=
0

,


a

3

n4


=


e

3

n




cos


(

α
+

γ
n

+

θ
n


)



cos





β


,


a

3

n5


=
0

,


a

3

n6


=

cos





α





cos





β







(
172
)













From the above, Lagrange's equation becomes















t




(



L





q
.

j



)


-



L




q
j




=


Q
j

+




l
,
n





λ
ln



a

ln





j









(
173
)













where






q


1


=β, q


2


=α, q


3i





i


, q


4i





i


, q


5i


=z


6i


, q


6i


=z


12i


q


3ii





ii


, q


4ii





ii


, q


5ii


=z


6ii


, q


6ii


=z


12ii


  (174)






can be obtained as follows:















t




(



L




β
.



)


-



L



β



=




F




β
.



+




l
,
n





λ
ln



a

ln





1









(
175
)






















β
¨





m
saw2n


+

m
baI

+



m
b



(



b
0


sin





α

+


c
0


cos





α


)


2

+


m
sn




{



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+


c

1

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}

2


+


m
an




{



e

1

n



sin






(

α
+

γ
n

+

θ
n


)


+


c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}

2


+


m
wn




{



e

3

n



sin






(

α
+

γ
n

+
θ

)


+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}

2










+




2



β
.






α
.




m
b



(



b
0


sin





α

+


c
0


cos





α


)








(



b
0


cos





α

-


c
0


sin





α


)


+


m
sn



{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}



{




z
.


6

n




cos


(

α
+

γ
n

+

η
n


)



-


(


α
.

+


η
.

n


)



z

6

n




sin


(

α
+

γ
n

+

η
n


)



-


α
.



[



c

1

n



sin






(

α
+

γ
n


)


-


b

2

n



cos





α


]



}


+


m
an



{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}



{



(


α
.

+


θ
.

n


)



e

1

n




cos


(

α
+

γ
n

+

θ
n


)



-


α
.



[



c

2

n




sin


(

α
+

γ
n


)



-


b

2

n



cos





α


]



}


+


m
wn



{



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}



{



(


α
.

+


θ
.

n


)



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



-


α
.



[



c

2

n




sin


(

α
+

γ
n


)



-


b

2

n



cos





α


]



}







-

α
¨





m
ba



(



b
0


cos





α

-


c
0


sin





α


)



+



α
.

2




m
ba



(



b
0


sin





α

+


c
0


cos





α


)



-



z
¨


6

n




m
sn







a

1

n




cos


(

α
+

γ
n

+

η
n


)



+




z
.


6

n




(


α
.

+


η
.

n


)








m
sn







a

1

n








sin
(





α
+

γ
n

+

η
n


)


+



η
¨

n



m
sn



z

6

n








a

1

n








sin
(





α
+

γ
n

+

η
n


)


+



η
.

n



m
sn




z
.


6

n




a

1

n




sin


(

α
+

γ
n

+

η
n


)



+




η
.

n



(


α
.

+


η
.

n


)




m
sn



z

6

n




a

1

n




cos


(

α
+

γ
n

+

η
n


)



-



θ
¨

n



m
aw1n



a

1

n








cos


(

α
+

γ
n

+

θ
n


)



+




θ
.

n



(


α
.

+


θ
.

n


)








m
aw1n







a

1

n








sin
(





α
+

γ
n

+

θ
n


)


+


α
¨



a

1

n




{



m
sawcn







sin


(

α
+

γ
n


)



-


m
sawbn


cos





α

+


m
sn



z

6

n




sin


(

α
+

γ
n

+

η
n


)



-


m
aw1n



cos


(

α
+

γ
n

+

θ
n


)




}


+


α
.



a

1

n




{



α
.



m
sawcn



cos


(

α
+

γ
n


)



+


α
.



m
sawbn


sin





α

+


(


α
.

+


η
.

n


)







m
sn







z

6

n








cos
(





α
+

γ
n

+

η
n


)


+


m
sn




z
.


6

n








sin
(





α
+

γ
n

+

η
n


)


+


(


α
.

+


θ
.

n


)







m
aw1n







sin
(





α
+

γ
n

+

θ
n


)



}


-

g


{



m
ba


cos





β

+



m
b



(



b
0


sin





α

+


c
0


cos





α


)



sin





β


}


-


n

i
,
ii







g


[



{



m
sn



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


m
aw1n



sin


(

α
+

γ
n

+

θ
n


)



+


m
sawcn



cos
(





α
+

γ
n


)


+


m
sawbn


sin





α


}






sin





β

+


m
sawan


cos





β


]






=


λ

3

n




[



-

{



z

12

n



cos





α

+


e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}



sin





β

+


a

1

n



cos





β


]









(
176
)






















β
¨



(


m
saw2n

+

m
baI

+


m
b



A
1
2


+


m
sn



B
1
2


+


m
an



B
2
2


+


m

2

n




B
3
2



)


+

2



β
.





[



α
.



m
b



A
1



A
2


+


m
sn



B
1



{




z
.


6

n




C

αγη





n



-


(


α
.

+


η
.

n


)



z

6

n




S

αγη





n



-


α
.



A
4



}


+


m
an



B
2



{



(


α
.

+


θ
.

n


)



e

1

n




C

αγθ





n



-


α
.



A
6



}


+


m
wn



B
3



{



(


α
.

+


θ
.

n


)



e

3

n




S

αγθ





n



-


α
.



A
6



}



]


-


α
¨



m
ba



A
2


+



α
.

2



m
ba



A
1


-



z
¨


6

n




m
sn



a

1

n




C

αγη





n



+

2




z
.


6

n




(


α
.

+


η
.

n


)




m
sn



a

1

n




S

αγη





n



+



η
¨

n



m
sn



z

6

n




a

1

n




S

αγη





n



+




η
.

n



(


2


α
.


+


η
.

n


)




m
sn



z

6

n




a

1

n




C

αγη





n



-



θ
¨

n



m
aw1n



a

1

n




C

αγ





θ





n



+




θ
.

n



(


2


α
.


+


θ
.

n


)




m
aw1n



a

1

n




S

αγθ





n



+


α
¨



a

1

n




{



m
sawcn



S

αγ





n



-


m
sawbn



C
α


+


m
sn



z

6

n




S

αγη





n



-


m
aw1n



C

αγ






θ
n





}


+



α
.

2



a

1

n




{



m
sawcn



C

αγ





n



+


m
sawbn



S
α


+


m
sn



z

6

n




C

αγ





η





n



+


m
aw1n



S

αγθ





n




}


-

g


[



m
ba



C
β


+


m
b



A
1



S
β


+


{



m
sn



z

6

n




C

αγη





n



+


m
aw1n



S

αγθ





n



+


m
sawcn



C

αγ





n



+


m
sawbn



S
α



}



S
β


+


m
sawan



C
β



]



=


λ

3

n




[



-

{



z

12

n




C
α


+


e

3

n




S

αγθ





n



+


c

2

n




C

αγ





n



+


b

2

n




S
α



}




S
β


+


a

1

n




C
β



]






(
177
)




















β
¨

=





2



β
.





[



α
.



m
b



A
1



A
2


+


m
sn



B
1



{




z
.


6

n




C

αγη





n



-


(


α
.

+


η
.

n


)



z

6

n




S

αγη





n



-














α
.



A
4


}

+


m
an



B
2



{



(


α
.

+


θ
.

n


)



e

1

n




C

αγθ





n



-


α
.



A
6



}


+









m
wn



B
3



{



(


α
.

+


θ
.

n


)



e

3

n




S

αγθ





n



-


α
.



A
6



}


]

-


α
¨



m
ba



A
2


+









α
.

2



m
ba



A
1


-



z
¨


6

n




m
sn



a

1

n




C

αγη





n



+

2




z
.


6

n




(


α
.

+


η
.

n


)




m
sn



a

1

n




S

αγη





n



+









η
¨

n



m
sn



z

6

n




a

1

n




S

αγη





n



+




η
.

n



(


2


α
.


+


η
.

n


)




m
sn



z

6

n




a

1

n




C

αγη





n



-









θ
¨

n



m
aw1n



a

1

n




C

αγθ





n



+




θ
.

n



(


2






α
.


+


θ
.

n


)




m
aw1n



a

1

n




S

αγθ





n



+







α
¨



a

1

n




{



m
sawcn



S

αγ





n



-


m
sawbn



C
α


+


m
sn



z

6

n




S

αγη





n



-











m
aw1n



C

αγθ





n



}

+



α
.

2



a

1

n




{



m
sawcn



C

αγ





n



+


m
sawbn



S
α


+













m
sn



z

6

n




C

αγη





n



+


m
aw1n



S

αγθ





n




}

-

g
[



m
ba



C
β


+










m
b



A
1



S
β


+

{



m
sn



z

6

n




C

αγη





n



+


m
aw1n



S

αγθ





n



+


m
sawcn



C
αγη


+














m
sawbn



S
α


}



S
β


+


m
sawan



C
β



]

+


λ

3

n




{

(



z

12

n




C
α


+
















e

3

n




S

αγθ





n



+


c

2

n




C

αγ





n



+


b

2

n




S
α



)



S
β


-


a

1

n




C
β



}





-

(


m
saw2n

+

m
baI

+


m
b



A
1
2


+


m
sn



B
1
2


+


m
an



B
2
2


+


m
wn



B
3
2



)







(
178
)












t




(



L








α
.



)


-



L



α



=




F




α
.



+




l
,
n





λ
ln



a

ln





2









(
179
)
























-

β
¨





m
ba



(



b
0


cos





α

-


c
0


sin





α


)



+


β
.



α
.




m
ba



(



b
0


sin





α

+


c
0


cos





α


)



+


α
¨







(


m
bbI

+

m
sawIn

+


m
sn




z

6

n




[


z

6

n


+

2


{



c

1

n



cos






η
n


-


b

2

n




sin


(


γ
n

+

η
n


)




}



]



-

2


m
aw1n



{



c

2

n



sin






θ
n


-


b

2

n



cos






(


γ
n

+

θ
n


)



}






+


α
.







m
sn





z
.


6

n




[


z

6

n


+

2


{



c

1

n



cos






η
n


-


b

2

n




sin


(


γ
n

+

η
n


)




}



]



+


m
sn




z

6

n




[



z
.


6

n


-

2



η
.

n



{



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+

η
n


)



}



]



-

2



θ
.

n



m
aw1n



{



c

2

n



cos






θ
n


+


b

2

n




sin


(


γ
n

+

θ
n


)




}






+



z
¨


6

n




m
sn



{



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+

η
n


)



}


+



z
.


6

n





η
.

n



m
sn



{



c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+



η
¨

n



m
sn



z

6

n




{


z

6

n


+


c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+



η
.

n



m
sn




z
.


6

n




{


z

6

n


+


c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+



η
.

n



m
sn



z

6

n




{



z
.


6

n


-



η
.

n



[



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+

η
n


)



]



}


+


θ
¨



[


m
aw2In

-


m
aw1n



{



c

2

n



sin






θ
n


-


b

2

n



cos






(


γ
n

+

θ
n


)



}



]


-



θ
.

n
2



m
aw1n



{



c

2

n



cos






θ
n


+


b

2

n



sin






(


γ
n

+

θ
n


)



}






]

+


β
¨



a

1

n




{



m
sawcn


sin






(

α
+

γ
n


)


-


m
sawbn


cos





α

+


m
sn



z

6

n



sin






(

α
+

γ
n

+

η
n


)


-


m
aw1n



cos


(

α
+

γ
n

+

θ
n


)




}


+


β
.



a

1

n




{



α
.



[



m
sawcn


cos






(

α
+

γ
n


)


+


m
sawbn


sin





α


]


+


m
sn




z
.


6

n



sin






(

α
+

γ
n

+

η
n


)


+


(


α
.

+


η
.

n


)



m
sn



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+


(


α
.

+


θ
.

n


)



m
aw1n



sin


(

α
+

γ
n

+

θ
n


)




}


-


{




β
.

2




m
b



(



b
0


cos





α

-


c
0


sin





α


)



+


α
.



β
.



m
ba



}



(



b
0


sin





α

+


c
0


cos





α


)


-



n

i
,
ii





&LeftBracketingBar;



+


β
.

2








m
sn



{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}



{



-

z

6

n





sin


(

α
+

γ
n

+

η
n


)



-


c

1

n



sin






(

α
+

γ
n


)


+


b

2

n



cos





α


}


+


m
an



{



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos
(





α
+

γ
n


)


+


b

2

n



sin





α


}



{



e

1

n




cos


(

α
+

γ
n

+

θ
n


)



-


c

2

n




sin


(

α
+

γ
n


)



+


b

2

n



cos





α


}


+


m
wn



{



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}



{



e
3



cos


(

α
+

γ
n

+

θ
n


)



-


c

2

n



sin






(

α
+

γ
n


)


+


b

2

n



cos





α


}






+



z
.


6

n




β
.



m
sn



a

1

n



sin






(

α
+

γ
n

+

η
n


)


+



η
.

n



β
.



m
sn







z

6

n








a

1

n








cos


(

α
+

γ
n

+

η
n


)



+


θ
.



β
.



m
aw1n







a

1

n








sin




(





α
+

γ
n

+

θ
n


)


+






α
.







β
.







a

1

n








{



m
sawcn



cos


(

α
+

γ
n


)



+


m
sawbn


sin





α

+


m
sn



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


m
aw1n


sin






(

α
+

γ
n

+

θ
n


)



}



&RightBracketingBar;




gm
b



(



b
0


cos





α

-


c
0


sin





α


)



cos





β


-







n

i
,
ii




g






{



m
sn



z

6

n




sin


(

α
+

γ
n

+

η
n


)



-


m
aw1n


cos






(

α
+

γ
n

+

θ
n


)


+


m
sawcn



sin


(

α
+

γ
n


)



-


m
sawbn


cos





α


}


cos





β



=


λ

3

n




{



-

z

12

n




sin





α

+


e

3

n



cos






(

α
+

γ
n

+

θ
n


)


-


c

2

n



sin






(

α
+

γ
n


)


+


b

2

n



cos





α


}


cos





β





(
180
)






















-

β
¨




m
ba



A
2


+


α
¨



{


m
bbI

+

m
saw1n

+


m
sn




z

6

n




(


z

6

n


+

2


E

1

n




)



-

2


m
aw1n



H

1

n




}


+

2


α
.



{



m
sn





z
.


6

n




(


z

6

n


+

E

1

n



)



-


m
sn



z

6

n





η
.

n



E

2

n



-



θ
.

n



m
aw1n



H

2

n




}


+



z
¨


6

n




m
sn



E

2

n



+



z
.


6

n





η
.

n



m
sn



E

1

n



+



η
¨

n



m
sn



z

6

n




{


z

6

n


+

E

1

n



}


+



η
.

n



m
sn




z
.


6

n




{


2


z

6

n



+

E

1

n



}


-



η
.

n
2



m
sn



z

6

n




E

2

n



+


θ
¨



(


m
aw2In

-


m
aw1n



H

1

n




)


-



θ
.

n
2



m
aw1n



H

2

n



+


β
¨




a

1

n




(



m
sawcn



S

αγ





n



-


m
sawbn



C
α


+


m
sn



z

6

n




S

αγη





n



-


m
aw1n



C

αγθ





n




)



+


β
.



a

1

n




{



α
.



(



m
sawcn



C

αγ





n



+


m
sawbn



S
α



)


+


m
sn




z
.


6

n




S

αγη





n



+


(


α
.

+


η
.

n


)







m
sn







z

6

n








C

αγη





n



+


(


α
.

+


θ
.

n


)



m
aw1n



S

αγθ





n




}


-



β
.

2



m
b



A
2



A
1


-

[



β
.



{



m
sn




B
1



(



-

z

6

n





S

αγη





n



-

A
4


)



+


m
an




B
2



(



e
1



C

αγθ





n



-

A
6


)



+


m
wn




B
3



(



e
3



C

αγθ





n



-

A
6


)




}


+



z
.


6

n




β
.



m
sn



a

1

n




S

αγη





n



+



η
.

n



β
.



m
sn



z

6

n




a

1

n




C

αγη





n



+


θ
.



β
.



m
aw1n



a

1

n




S

αγθ





n



+


α
.



β
.



a

1

n




{



m
sawcn



C

αγ





n



+


m
sawbn



S
α


+


m
sn



z

6

n




C

αγ





η





n



+


m
aw1n



S

αγθ





n




}



]

+


gm
b



A
2



C
β


-

g


{



m
sn



z

6

n




S

αγη





n



-


m
aw1n



C

αγθ





n



+


m
sawcn



S

αγ





n



-


m
sawbn



C

α








}



C
β



=


λ

3

n




{



-

z

12

n





S
α


+


e

3

n




C

αγ





θ





n



-


c

2

n




S

αγ





n



+


b

2

n




C
α



}



C
β






(
181
)






















-

β
¨




m
ba



A
2


+


α
¨



{


m
bbI

+

m
saw1n

+


m
sn




z

6

n




(


z

6

n


+

2


E

1

n




)



-

2


m
aw1n



H

1

n




}


+



m
sn



(


2


α
.




z
.


6

n



+



η
¨

n



z

6

n



+

2


η
.




z
.


6

n




)




(


z

6

n


+

E

1

n



)


-

2



α
.



(



m
sn



z

6

n





η
.

n



E

2

n



+



θ
.

n



m
aw1n



H

2

n




)



+



z
¨


6

n




m
sn



E

2

n



-



η
.

n
2



m
sn



z

6

n




E

2

n



+


θ
¨



(


m
aw2In

-


m
aw1n



H

1

n




)


-



θ
.

n
2



m
aw1n



H

2

n



+


β
¨




a

1

n




(



m
sawcn



S

αγ





n



-


m
sawbn



C
α


+


m
sn



z

6

n




S

αγη





n



-


m
aw1n



C

αγθ





n




)



-



β
.

2



{



m
b



A
2



A
1


+


m
sn




B
1



(



-

z

6

n





S

αγη





n



-

A
4


)



+


m
an




B
2



(



e
1



C

αγθ





n



-

A
6


)



+


m
wn




B
3



(



e
3



C

αγθ





n



-

A
6


)




}


+


gm
b



A
2



C
β


-

g


{



m
sn



z

6

n




S

αγη





n



-


m
aw1n



C

αγθ





n



+


m
sawcn



S

αγ





n



-


m
sawbn



C
α



}



C
β



=



λ

3

n




(



-

z

12

n





S
α


+


e

3

n




C

αγθ





n



-


c

2

n




S

αγ





n



+


b

2

n




C
α



)




C
β






(
182
)

















α
¨

=








m
sn



(


2


α
.




z
.


6

n



+



η
¨

n



z

6

n



+

2



η
.

n




z
.


6

n




)




(


z

6

n


+

E

1

n



)


-







2



α
.



(



m
sn



z

6

n





η
.

n



E

2

n



+



θ
.

n



m
aw1n



H

2

n




)



+



z
¨


6

n




m

sn








E

2

n



-









η
.

n
2



m
sn



z

6

n




E

2

n



+


θ
¨



(


m
aw2In

-


m
aw1n



H

1

n




)


-









θ
.

n
2



m
aw1n



H

2

n



+


β
¨



a

1

n


(



m
sawcn



S

αγ





n



-


m
sawbn



C
α


+












m
sn



z

6

n




S

αγη





n



-


m
aw1n



C

αγθ





n




)

-



β
.

2



{



m
b



A
2



A
1


+











m
sn




B
1



(



-

z

6

n





S

αγη





n



-

A
4


)



+


m
an




B
2



(



e
1



C

αγθ





n



-

A
6


)



+









m
wn




B
3



(



e
3



C

αγθ





n



-

A
6


)



}

+


gm
b



A
2



C
β


-






g


{



m
sn



z

6

n




S

αγη





n



-


m
aw1n



C

αγθ





n



+


m
sawcn



S

αγ





n



-












m
sawbn



C
α


}



C
β


-


β
¨



m
ba



A
2


+


λ

3

n


(



z

12

n




S
α


-












e

3

n




C

αγθ





n



+


c

2

n




S

αγ





n



-


b

2

n




C
α



)



C
β






-

{


m
bbI

+

m
sawIn

+


m
sn




z

6

n




(


z

6

n


+

2


E

1

n




)



-

2


m
aw1n



H

1

n




}







(
183
)


























t




(



L





η
.

n



)


-



L




η
n




=




F





η
.

n



+




l
,
n





λ
ln



a

ln





3













l
=
1

,
2
,


3





n

=
i

,
ii





(
184
)









m
sn




η
¨

n



z

6

n

2


+

2


m
sn




η
.

n




z
.


6

n




z

6

n



+


α
¨



m
sn



z

6

n




{


z

6

n


+


c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+


α
.



m
sn




z
.


6

n




{


z

6

n


+


c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+


α
.



m
sn



z

6

n




{



z
.


6

n


-



η
.

n



[



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+

η
n


)



]



}


+


β
¨



m
sn



z

6

n




a

1

n



sin






(

α
+

γ
n

+

η
n


)


+


β
.



m
sn




z
.


6

n




a

1

n



sin






(

α
+

γ
n

+

η
n


)


+



β
.



(


α
.

+


η
.

n


)




m
sn



z

6

n




a

1

n



cos






(

α
+

γ
n

+

η
n


)


-






α
.

2



m
sn



z

6

n




{



-

c

1

n




sin






η
n


-


b

2

n



cos






(


γ
n

+

η
n


)



}


+



β
.

2



m
sn



{



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+


c

1

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}



{


-

z

6

n




sin






(

α
+

γ
n

+

η
n


)


}


+



z
.


6

n




α
.



m
sn



{



c

1

n



cos






η
n


-


b

2

n



sin






(


γ
n

+

η
n


)



}


+



z
.


6

n




β
.



m
sn



a

1

n




sin




(

α
+

γ
n

+

η
n


)


-



η
.

n



α
.



m
sn



z

6

n




{



c

1

n



sin






η
n


+


b

2

n



cos






(


γ
n

+

η
n


)



}


+



η
.

n



β
.



m
sn







z

6

n








a

1

n








cos




(

α
+

γ
n

+

η
n


)


+


α
.



β
.



a

1

n




m
sn



z

6

n



cos






(

α
+

γ
n

+

η
n


)


+


gm
sn



z

6

n




sin


(

α
+

γ
n

+

η
n


)



cos





β





=



-


λ

1

n




(


z

6

n


-

d

1

n



)




cos






η
n


+



λ

2

n




(


z

6

n


-

d

1

n



)



sin






η
n







(
185
)






















m
sn




η
¨

n



z

6

n

2


+

2


m
sn




η
.

n




z
.


6

n



+


α
¨



m
sn



z

6

n




{


z

6

n


+

E
1


}


+


α
.



m
sn




z
.


6

n




{


2


z

6

n



+

E
1


}


-


α
.



m
sn



z

6

n





η
.

n



E
2


+


β
¨



m
sn



z

6

n




a

1

n




S

αγη





n



+


β
.



m
sn




z
.


6

n




a

1

n




S

αγη





n



+



β
.



(


α
.

+


η
.

n


)




m
sn



z

6

n




a

1

n




C

αγη





n



+



α
.

2



m
sn



z

6

n




E
2


+



β
.

2



m
sn



B
1



z

6

n




S

αγη





n



-



z
.


6

n




α
.



m
sn



E
1


-



z
.


6

n




β
.



m
sn



a

1

n




S

αγη





n



+



η
.

n



α
.



m
sn



z

6

n




E
2


-



η
.

n



β
.



m
sn



z

6

n




a

1

n




C

αγη





n



-


α
.



β
.



a

1

n




m
sn



z

6

n




C

αγη





n



-


gm
sn



z

6

n




S

αγη





n




C
β



=



-


λ

1

n




(


z

6

n


-

d

1

n



)





C

η





n



+



λ

2

n




(


z

6

n


-

d

1

n



)




S

η
n








(
186
)








m
sn



z

6

n




{




η
¨

n



z

6

n



+

2



η
.

n




z
.


6

n



+


α
¨



(


z

6

n


+

E
1


)


+

2


α
.




z
.


6

n



+


β
¨



a

1

n




S

αγη





n



+



α
.

2



E
2


+



β
.

2



B
1



S

αγη





n



-


gS

αγη





n




C
β



}


=



-


λ

1

n




(


z

6

n


-

d

1

n



)





C

η





n



+



λ

2

n




(


z

6

n


-

d

1

n



)




S

η





n








(
187
)

















λ

1

n


=






m
sn



z

6

n




{




η
¨

n



z

6

n



+

2



η
.

n




z
.


6

n



+


α
¨



(


z

6

n


+

E
1


)


+











2






α
.




z
.


6

n



+


β
¨



a

1

n




S

αγη





n



+



α
.

2



E
2


+



β
.

2



B
1



S

αγη





n



-


gS

αγη





n




C
β



}

-








λ

2

n




(


z

6

n


-

d

1

n



)




S

η





n








-

(


z

6

n


-

d

1

n



)




C

η





n








(
188
)


























t




(



L





θ
.

n



)


-



L




θ
n




=




F





θ
.

n



+




l
,
n





λ

1

n




a

ln





4

















l
=
1

,
2
,


3





n

=
i

,
ii





(
189
)









θ
¨

n



m
aw2In


+


α
¨

[




m
aw2In

-


m
aw1n



{



c

2

n



sin






θ
n


-


b

2

n



cos






(


γ
n

+

θ
n


)



}


-


α
.




θ
.

n



m
aw1n



{



c

2

n



cos






θ
n


+


b

2

n



sin






(


γ
n

+

θ
n


)



}


-


β
¨



m
aw1n



a

1

n




cos


(

α
+

γ
n

+

θ
n


)



+



β
.



(


α
.

+


θ
.

n


)




m
aw1n



a

1

n








sin




(

α
+

γ
n

+

θ
n


)






-



[



k
zi



e

0

i

2



{


sin


(


γ
i

+

θ
i


)


-

sin


(


γ
ii

+

θ
ii


)



}



cos


(


γ
n

+

θ
n


)




X
s


-



α
.

2



m
aw1n



{



c

2

n



cos






θ
n


+


b

2

n




sin


(


γ
n

+

θ
n


)




}


+


β
.

2











m
an







{



e

1

n








sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}







e

1

n




cos


(

α
+

γ
n

+

θ
n


)



+


m
wn



{



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}







e

3

n








cos
(





α
+

γ
n

+

θ
n


)







-

θ
.




α
.



m
aw1n



{



c

2

n



cos






θ
n


+


b

2

n



sin






(

γ
+

θ
n


)



}


+


θ
.



β
.



m
aw1n







a

1

n








sin
(





α
+

γ
n

+

θ
n


)


+


α
.



β
.



a

1

n




m
aw1n



sin


(

α
+

γ
n

+

θ
n


)



-


gm
aw1n



cos


(

α
+

γ
n

+

θ
n


)



cos





β





]

=



λ

1

n




e

2

n



sin






θ
n


+


λ

2

n




e

2

n



cos






θ
n


+


λ

3

n




e

3

n



cos






(

α
+

γ
n

+

θ
n


)


cos





β








(
190
)























θ
¨

n



m
aw2In


+


α
¨



(


m
aw2In

-


m
aw1n



H
1



)


-


α
.




θ
.

n



m
aw1n



H
2


-


β
¨



m
aw1n



a

1

n




C

αγθ





n



+



β
.



(


α
.

+


θ
.

n


)








m
aw1n







a

1

n








S

αγθ





n



-

[



-

k
zi








e

0

i

2







{


sin


(


γ
i

+

θ
i


)


+

sin


(


γ
ii

+

θ
ii


)



}







cos




(






γ
n

+

θ
n


)



X
s






-



α
.

2



m
aw1n



H
2


+



β
.

2



(



m
an



B
2



e

1

n




C

αγθ





n



+


m
wn



B
3



e

3

n




C

αγθ





n




)


-


θ
.



α
.



m
aw1n



H
2


+


θ
.



β
.



m
aw1n



a

1

n




S

αγθ





n



+


α
.



β
.



a

1

n




m
aw1n



S

αγθ





n



-


gm
aw1n



C

αγθ





n




C
β



]


=



λ

1

n




e

2

n




S

θ





n



+


λ

2

n




e

2

n




C

θ





n



+


λ

3

n




e

3

n




C

αγθ





n




C
β







(
191
)










θ
¨

n



m
aw2In


+


α
¨



(


m
aw2In

-


m
aw1n



H
1



)


-


β
¨



m
aw1n



a

1

n




C

αγθ





n



+



α
.

2



m
aw1n



H
2


-



β
.

2



(



m
an



B
2



e

1

n




C

αγθ





n



+


m
wn



B
3



e

3

n




C

αγθ





n




)


+


gm
aw1n



C

αγθ





n




C
β


+


k
zi



e

0

i

2



{


sin


(


γ
i

+

θ
i


)


+

sin


(


γ
ii

+

θ
ii


)



}


cos






(


γ
n

+

θ
n


)



=



λ

1

n




e

2

n




S

θ





n



+


λ

2

n




e

2

n




C

θ





n



+


λ

3

n




e

3

n




C

αγθ





n




C
β







(
192
)































θ
¨

n

=







α
¨



(


m
aw2In

-


m
aw1n



H
1



)


-


β
¨



m
aw1n



a

1

n




C

αγθ





n



+









α
.

2



m
aw1n



H
2


-



β
.

2

(



m
an



B
2



e

1

n




C

αγθ





n



+











m
wn



B
3



e

3

n




C

αγθ





n



)

+


gm
aw1n



C

αγθ





n




C
β


-








λ

1

n




e

2

n




S

θ





n



-


λ

2

n




e

2

n




C

θ





n



-


λ

3

n




e

3

n




C

αγθ





n




C
β


+







k
zi



e

0

i

2



{


sin


(


γ
i

+

θ
i


)


+

sin


(


γ
ii

+

θ
ii


)



}


cos






(


γ
n

+

θ
n


)






-

m
aw2In







(
193
)













t




(



L





z
.


6

n




)


-



L




z

6

n





=




F





z
.


6

n




+




l
,
n





λ
ln



a

ln





5













l
=
1

,
2
,


3





n

=
i

,
ii





(
194
)






















m
sn




z
¨


6

n



+


α
¨



m
sn



{



c

1

n



sin






η
n


+


b

2

n




cos


(


γ
n

+

η
n


)




}


+


α
.




η
.

n



m
sn



{



c

1

n



cos






η
n


-


b

2

n




sin


(


γ
n

+

η
n


)




}


-


β
¨



m
sn



a

1

n




cos


(

α
+

γ
n

+

η
n


)



+



β
.



(


α
.

+


η
.

n


)




m
sn



a

1

n




sin


(

α
+

γ
n

+

η
n


)



-





m
sn




η
.

n
2



z

6

n



+



α
.

2




m
sn



[


z

6

n


+

{



c

1

n



cos






η
n


-


b

2

n




sin


(


γ
n

+

η
n


)




}


]



+



β
.

2



m
sn



{



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}







cos
(





α
+

γ
n

+

η
n


)


+



η
.

n



α
.



m
sn







{


2


z

6

n



+


c

1

n



cos






η
n


-


b

2

n




sin


(


γ
n

+

η
n


)




}


+



η
.

n







β
.







m
sn







a

1

n








sin


(

α
+

γ
n

+

η
n


)



+


α
.



β
.



a

1

n




m
sn



sin


(

α
+

γ
n

+

η
n


)



-


gm
sn



cos


(

α
+

γ
n

+

η
n


)



cos





β

-


k
sn



(


z

6

n


-

l
sn


)






=



-

c
sn





z
.


6

n



-


λ

1

n



sin






η
n


-


λ

2

n



cos






η
n







(
195
)









m
sn



{



z
¨


6

n


+


α
¨



E
2


-


β
¨



a

1

n




C

αγη





n



-



η
.

n
2



z

6

n



-



α
.

2



(


z

6

n


+

E
1


)


-



β
.

2



B
1



C

αγθ





n



-

2



η
.

n



α
.



z

6

n



+


gC

αγη





n




C
β



}


+


k
sn



(


z

6

n


-

l
sn


)



=



-

c
sn





z
.


6

n



-


λ

1

n




S

η





n



-


λ

2

n




C

η





n








(
196
)


















λ

2

n


=






m
sn



{



z
¨


6

n


+


α
¨



E
2


-


β
¨



a

1

n




C

αγη





n



-



η
.

n
2



z

6

n



-



α
.

2



(


z

6

n


+

E
1


)


-













β
.

2



B
1



C

αγη





n



-

2



η
.

n



α
.



z

6

n



+


gC

αγη





n




C
β



}

+








k
sn



(


z

6

n


-

l
sn


)


+


c
sn




z
.


6

n



+


λ

1

n




S

η





n








-

C

η





n





















t




(



L





z
.


12

n




)


-



L




z

12

n





=




F





z
.


12

n




+




l
,
n





λ
ln



a

ln





6













l
=
1

,
2
,


3





n

=
i

,
ii













k
wn



(


z

12

n


-

l
wn


)


=



-

c
wn





z
.


12

n



+


λ

3

n



cos





αcos





β








=



-

c
wn





z
.


12

n



+


λ

3

n




C
α



C
β

















(
198
)







λ

3

n


=




c

2

n





z
.


12

n



+


k
wn



(


z

12

n


-

l
wn


)




C
α






(
199
)








(
197
)













From the differentiated constraints it follows that:






{dot over (θ)}


n


e


2n


sin θ


n


−{dot over (z)}


6n


sin {dot over (η)}


n


(z


6n


−d


ln


)cos ηn=0








{dot over (θ)}


n


e


2n


cos θ


n


−{dot over (z)}


6n


cos η


n


+{dot over (η)}


n


(z


6n


−d


1n


)sin η


n


=0  (200)






and












{



z

12

n



cos





α

-


α
.



z

12

n



sin





α

+


(


α
.

+


θ
.

n


)



e

3

n




cos


(

α
+

γ
n

+

θ
n


)



-


α
.



c

2

n




sin


(

α
+

γ
n


)



+


α
.



b

2

n



cos





α


}


cos





β

-


β
.



[



{



z

12

n



cos





α

+


e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


sin





β

+


a

1

n



cos





β


]


-



R
.

n



(
t
)



=
0




(
201
)













Second order of differentiation of equation (200) yields:














θ
¨

n



e

2

n




S

θ





n



+



θ
n
2

.



e

2

n




C

θ





n



-



z
¨


6

n




S

η





n



-



z
.


6

n





η
.

n



C

η





n



-




η
¨

n



(


z

6

d


-

d

1

n



)




C

η





n



-



η
.

n




z
.


6

n




C

η





n



+




η
.

n
2



(


z

6

n


-

d

1

n



)




S

η





n




=
0











θ
¨

n



e

2

n




C

θ





n



+



θ
n
2

.



e

2

n




S

θ





n



-



z
¨


6

n




C

η





n



-



z
.


6

n





η
.

n



S

η





n



-




η
¨

n



(


z

6

d


-

d

1

n



)




S

η





n



-



η
.

n




z
.


6

n




S

η





n



+




η
.

n
2



(


z

6

n


-

d

1

n



)




C

η





n




=
0





(
202
)


















η
¨

n

=





θ
¨

n



e

2

n




S

θ





n



+



θ
.

n
2



e

2

n




C

θ





n



-



z
¨


6

n




S

η





n



-

2



η
.

n




z
.


6

n




C

η





n



+




η
.

n
2



(


z

6

n


-

d

1

n



)




S

η





n






(


z

6

n


-

d

1

n



)



C

η





n








(
203
)









z
¨


6

n


=





θ
¨

n



e

2

n




C

θ





n



+



θ
.

n
2



e

2

n




S

θ





n



+




η
¨

n



(


z

6

n


-

d

1

n



)




S

η





n



+

2



η
.

n




z
.


6

n




S

η





n



+




η
.

n
2



(


z

6

n


-

d

1

n



)




C

η





n





C

η





n









and




(
204
)








z
.


12





n


=







{



α
.



z

12

n




S
α


-


(


α
.

+


θ
.

n


)



e

3

n




C

αγ





θ





n



+


α
.



c

2

n




S

αγ





n



-


α
.



b

2

n




C
α



}



C
β


+








β
.



[



{



z

12

n




C
α


+


e

3

n




S

αγ





θ





n



+


c

2

n




C

αγ





n



+


b

2

n




S
α



}



S
β


+


a

1

n




C
β



]






R
.

n



(
t
)








C
α



C
β







(
205
)













Supplemental differentiation of equation (199) for the later entropy production calculation yields:






k


wn


{dot over (z)}


12n


=−c


wn


{umlaut over (z)}


12n


+{dot over (λ)}


3n


C


α


C


β


−{dot over (α)}λ


3n


S


α


C


β


−{dot over (β)}λ


3n


C


α


C


β


  (206)






therefore











z
¨


12

n


=





λ
.


3

n




C
α



C
β


-


α
.



λ

3

n




S
α



C
β


-


β
.



λ

3

n




C
α



S
β


-


k

2

n





z
.


12

n





c
wn






(
207
)













or from the third equation of constraint:












{




z
¨


12

n



cos





α

-



z
.


12

n




α
.


cos





α

-


α
¨



z

12

n



sin





α

-


α
.




z
.


12

n



sin





α

-



α
.

2



z

12

n



cos





α

+


(


α
¨

+


θ
.

n


)



e

3

n



cos






(

α
+

γ
n

+

θ
n


)


-



(


α
.

+


θ
.

n


)

2



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



-


α
¨



c

2

n



sin






(

α
+

γ
n


)


-



α
.

2



c

2

n



cos






(

α
+

γ
n


)


+


α
¨



b

2

n



cos





α

-



α
.

2



b

2

n



sin





α


}


cos





β

-


β
.



{




z
.


12

n



cos





α

-


α
.



z

12

n



sin





α

+


(


α
.

+


θ
.

n


)



e

3

n



cos






(

α
+

γ
n

+

θ
n


)


-


α
.



c

2

n




sin


(

α
+

γ
n


)



+


α
.



b

2

n



cos





α


}


sin





β

-


β
¨



[



{



z

12

n



cos





α

+


e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}






sin





β

+


a

1

n



cos





β


]


-


β
.



[



{




z
.


12

n



cos





α

-


α
.



z

12

n



sin





α

+


(


α
.

+


θ
.

n


)



e

3

n




cos


(

α
+

γ
n

+

θ
n


)



-


(


α
.

+


γ
.

n


)



c

2

n




sin


(

α
+

γ
n


)



+


α
.



b

2

n



cos





α


}






sin





β

+


β
.



{



z

12

n



cos





α

+


e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}






cos





β

-


β
.



a

1

n



sin





β


]


-



R
¨

n



(
t
)



=
0




(
208
)








z
¨


12

n


=





{



-


z
.


12

n





α
.


cos





α

-


α
¨



z

12

n



sin





α

-


α
.




z
.


12

n



sin





α

-



α
.

2



z

12

n



cos





α

+









(


α
¨

+


θ
.

n


)



e

3

n



cos






(

α
+

γ
n

+

θ
n


)


-



(


α
.

+


θ
.

n


)

2



e

3

n



sin






(

α
+

γ
n

+

θ
n


)


-








α
¨



c

2

n



sin






(

α
+

γ
n


)


-



α
.

2



c

2

n



cos






(

α
+

γ
n


)


+


α
¨



b

2

n



cos





α

-











α
.

2



b

2

n



sin





α

}


cos





β

-


β
.



{




z
.


12

n



cos





α

-


α
.



z

12

n



sin





α

+











(


α
.

+


θ
.

n


)



e

3

n



cos






(

α
+

γ
n

+

θ
n


)


-


α
.



c

2

n



sin






(

α
+

γ
n


)


+










α
.



b

2

n



cos





α

}


sin





β

-


β
¨

[

{



z

12

n



cos





α

+


e

3

n



sin






(

α
+

γ
n

+

θ
n


)


+
















c

2

n



cos






(

α
+

γ
n


)


+


b

2

n



sin





α


}


sin





β

+


a

1

n



cos





β


]

-







β
.

[

{




z
.


12

n



cos





α

-


α
.



z

12

n



sin





α

+


(


α
.

+


θ
.

n


)



e

3

n



cos






(

α
+

















γ
n

+

θ
n


)

-


(


α
.

+


γ
.

n


)



c

2

n




sin


(

α
+

γ
n


)



+


α
.



b

2

n



cos





α


}


sin





β

+







β
.



{



z

12

n



cos





α

+


e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+













c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α


}


cos





β

-









β
.



a

1

n



sin





β

]

-



R
¨

n



(
t
)








-
cos






αcos





β






(
209
)













5. Summarization for simulation.










β
¨

=





2


β
.





[



α
.



m
b



A
1



A
2


+


m
sn



B
1



{




z
.


6

n




C

αγ





η





n



-


(


α
.

+


η
.

n


)



z

6

n




S

αγη





n



-













α
.



A
4


}

+


m
an



B
2



{



(


α
.

+


θ
.

n


)



e

1

n




C

αγθ





n



-


α
.



A
6



}


+









m
wn



B
3



{



(


α
.

+


θ
.

n


)



e

3

n




S

αγθ





n



-


α
.



A
6



}


]

-








α
¨



m
ba



A
2


+



α
.

2



m
ba



A
1


-



z
¨


6

n




m
sn



a

1

n




C

αγ





η





n



+







2




z
.


6

n




(


α
.

+


η
.

n


)




m
sn



a

1

n




S

αγη





n



+



η
¨

n



m
sn



z

6

n




a

1

n




S

αγη





n



+










η
.

n



(


2


α
.


+


η
.

n


)




m
sn



z

6

n




a

1

n




C

αγ





η





n



-









θ
¨

n



m
aw1n



a

1

n




C

αγθ





n



-




θ
.

n



(


2


α
.


+


θ
.

n


)




m
aw1n



a

1

n




S

αγθ





n



+







α
¨



a

1

n




{



m
sawcn



S

αγ





n



-


m
sawbn



C
α


+


m
sn



z

6

n




S

αγ





η





n



-











m
aw1n



C

αγθ





n



}

+



α
.

2



a

1

n




{



m
sawcn



C

αγ





n



+


m
sawbn



S
α


+













m
sn



z

6

n




C

αγη





n



+


m
aw1n



S

αγ





θ





n




}

-

g
[



m
ba



C
β


+










m
b



A
1



S
β


+

{



m
sn



z

6

n




C

αγη





n



+


m
aw1n



S

αγθ





n



+















m
sawcn



C

αγ





n



+


m
sawbn



S
α



}



S
β


+


m
sawan



C
β



]

-







λ

3

n




{

(



z

12

n




C
α


+


e

3

n




S

αγ





θ





n



+


c

2

n




C

α





γ





n



+














b

2

n




S
α


)



S
β


-


a

1

n




C
β



}






-

(


m
saw2n

+

m
ba1

+


m
b



A
1
2


+


m
sn



B
1
2


+


m
an



B
2
2


+


m

2

n




B
3
2



)






*
Sign





of






λ

3

n







is





modified





judging





from





the





behavior





of





the






model
.







(
210
)







α
¨

=








m
sn



(


2


α
.




z
.


6

n



+



η
¨

n



z

6

n



+

2



η
.

n




z
.


6

n




)








(


z

6

n


+

E

1

n



)


-







2



α
.



(



m
sn



z

6

n





η
.

n



E

2

n



+



θ
.

n



m
aw1n



H

2

n




)



+









z
¨


6

n




m
sn



E

2

n



-



η
.

n
2



m
sn



z

6

n




E

2

n



+


θ
¨

(


m
aw2In

-











m
aw1n



H

1

n



)

-



θ
n
2

.



m
aw1n



H

2

n



+


β
¨




a

1

n


(



m
sawcn



S

αγ





n



-













m
sawbn



C
α


+


m
sn



z

6

n




S

αγ





η





n



-


m
aw1n



C

αγ





θ





n




)

-








β
.

2



{



m
b



A
2



A
1


-


m
sn




B
1



(



z

6

n




S

αγη





n



+

A
4


)



+










m
an




B
2



(



e

1

n




C

αγθ





n



-

A
6


)



+


m
wn




B

3

n


(



e
3



C

αγθ





n



-












A
6

)

}

+


gm
b



A
2



C
β


-

g


{



m
sn



z

6

n




S

αγη





n



-


m
aw1n



C

αγθ





n



+














m
sawcn



S

αγ





n



-


m
sawbn



C
α



}



C
β


-


β
¨



m
ba



A
2


+








λ

3

n




(



z

12

n




S
α


-


e

3

n




C

αγθ





n



+


c

2

n




S

αγ





n



-


b

2

n




C
α



)




C
β






-

{


m
bbI

+

m
sawIn

+


m
sn




z

6

n




(


z

6

n


+

3


E

1

n




)



-

2


m
aw1n



H

1

n




}







(
211
)







λ

1

n


=






m
sn



z

6

n




{




η
¨

n



z

6

n



+

2



η
.

n




z
.


6

n



+


α
¨



(


z

6

n


+

E
1


)


+

2


α
.




z
.


6

n



+












β
¨



a

1

n




S

αγ





η





n



+



α
.

2



E
2


+



β
.

2



B
1



S

αγη





n



-


gS

αγ





η





n




C
β



}

-








λ

2

n




(


z

6

n


-

d

1

n



)




S

η





n








-

(


z

6

n


-

d

1

n



)




C

η





n








(
212
)








θ
¨

n

=







α
¨



(


m
aw2In

-


m
aw1n



H
1



)


-


β
¨



m
aw1n



a

1

n




C

αγθ





n



+



α
.

2



m
aw1n



H
2


-









β
.

2



(



m
an



B
2



e

1

n




C

αγθ





n



+


m
wn



B
3



e

3

n




C

αγθ





n




)


+


gm
aw1n



C

αγθ





n




C
β


-








λ

1

n




e

2

n




S


θ





n








-


λ

2

n




e

2

n




C

θ





n



-


λ

3

n




e

3

n




C

αγθ





n




C
β


+







1
2



k
zi



e

0

i

2



{


sin






(


γ
i

+

θ
i


)


+

sin






(


γ
ii

+

θ
ii


)



}


cos






(


γ
n

+

θ
n


)






-

m
aw2In







(
213
)













In the above expression, the potential energy of the stabilizer is halved because there are left and right portions of the stabilizer.










λ

2

n


=







m
sn





{
z

¨


6

n



+


α
¨



E
2


-


β
¨



a

1

n




C

αgη





n



-



η
.

n
2



z

6

n



-









α
.

2



(


z

6

n


+

E
1


)


-



β
.

2



B
1



C

αγη





n



-

2



η
.

n



α
.



z

6

n



+









gC

αγη





n




C
β


}

+


k
sn



(


z

6

n


+

l
sn


)


+


c

2

n





z
.


6

n



+


λ

1

n




S

η





n








-

C

η





n








(
214
)







λ

3

n


=




c
wn




z
.


12

n



+


k
wn



(


z

12

n


-

l
wn


)




C
α






(
215
)








η
¨

n

=








θ
¨

n



e

2

n




S

θ





n



+



θ
.

n
2



e

2

n




C

θ





n



-



z
¨


6

n




S

η





n



-

2



η
.

n




z
.


6

n




C

η





n



+









η
.

n
2



(


z

6

n


-

d

1

n



)




S

η





n








(


z

6

n


-

d

1

n



)



C

η





n








(
216
)








z
¨


6

n


=








θ
¨

n



e

2

n




C

θ





n



-



θ
.

n
2



e

2

n




S

θ





n



+




η
¨

n



(


z

6

n


-

d

1

n



)




S

η





n



+







2



η
.

n




z
.


6

n




S

η





n



+




η
.

n
2



(


z

6

n


-

d

1

n



)




C

η





n








C

η





n







(
217
)








z
.


12

n


=





{



α
.



z

12

n




S
α


-


(


α
.

+


θ
.

n


)



e

3

n




C

αγθ





n



+


α
.



c

2

n




S

α





γ





n



-











α
.



b

2

n




C
α


}



C
β


+


β
.

[

{



z

12

n




C
α


+


e

3

n




S

α





γθ





n



+
















c

2

n




C

αγ





n



+


b

2

n




S
α



}



S
β


+


a

1

n




C
β



]

+



R
.

n



(
t
)








C
α



C
β







(
218
)













where























































n
=
i

,
ii







m
ba

=


m
b



(


a
0

+

a
1


)












m
bbI

=



m
b



(


b
0
2

+

c
0
2


)


+

I
bx












m
baI

=




m
b



(


a
0

+

a
1


)


2

+

I
by












m
sawan

=


(


m
sn

+

m
an

+

m
wn


)



a

1

n













m
sawbn

=


(


m
sn

+

m
an

+

m
wn


)



b

2

n













m
sawcn

=



m
sn



c

1

n



+


(


m
an

+

m
wn


)



c

2

n














m
saw2n

=


(


m
sn

+

m
an

+

m
wn


)



a

1

n

2












m
sawIn

=



m
an



e

1

n

2


+


m
wn



e

3

n

2


+


m
sn



(


c

1

n

2

+

b

2

n

2

-

2


c

1

n




b

2

n



sin






γ
n



)


+












(


m
an

+

m
wn


)



(


c

2

n

2

+

b

2

n

2

-

2


c

2

n




b

2

n



sin






γ
n



)


+

I
axn











m
aw2In

=



m
an



e

1

n

2


+


m
wn



e

3

n

2


+

I
axn












m
aw1n

=



m
an



e

1

n



+


m
wn



e

3

n














m
aw2n

=



m
an



e

1

n

2


+


m
wn



e

3

n

2













A
1

=



b
0


sin





α

+


c
0


cos





α









A
2

=



b
0


cos





α

-


c
0


sin





α









A

4

n


=



c

1

n




sin


(

α
+

γ
n


)



-


b

2

n



cos





α









A

6

n


=



c

2

n




sin


(

α
+

γ
n


)



-


b

2

n



cos





α












B

1

n


=



z

6

n




cos


(

α
+

γ
n

+

η
n


)



+


c

1

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α












B

2

n


=



e

1

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α









B

3

n


=



e

3

n




sin


(

α
+

γ
n

+

θ
n


)



+


c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin





α









E

1

n


=



c

1

n



cos






η
n


-


b

2

n




sin


(


γ
n

+

η
n


)











E

2

n


=



c

1

n



sin






η
n


+


b

2

n




cos


(


γ
n

+

η
n


)











H

1

n


=



c

2

n



sin






θ
n


-


b

2

n




cos


(


γ
n

+

θ
n


)











H

2

n


=



c

2

n



cos






θ
n


+


b

2

n




sin


(


γ
n

+

θ
n


)











(
219
)









 S


α


=sin α, S


β


=sin β, S


αγn


=sin(α+γ


n


),






S


αγηn


=sin(α+γ


n





n


), S


αγθn


=sin(α+γ


n





n


)








C


α


=cos α, C


β


=cos β, C


αγn


=cos(α+γ


n


),








C


αγηn


=cos(α+γ


n





n


),








C


αγθ


=cos(α+γ


n





n


)  (220)






The initial conditions are:










β
.

=


α
.

=



η
.

n

=



θ
.

n

=



z
.


6

n


=



z
.


12

n


=
0









(
221
)






β
=

α
=
0





(
222
)











z

6

n


=






l
sn

+

l
s0n














l
s0n





0.6






m
b



2


k
sn




g













z

12

n


=






l
wn

+

l
w0n


















l
w0n





0.6


(


m
b

+

m
sn

+

m
an

+

m
wn


)



2


k
wn




g






(
223
)







η
n

=


cos

-
1







(


z

6

n


-

d

1

n



)

2

+


(


c

1

n


-

c

2

n



)

2

-

e

2

n

2




-
2



(


c

1

n


-

c

2

n



)







(


z

6

n


-

d

1

n



)








(
224
)







θ
n

=


sin

-
1







(


z

6

n


-

d

1

n



)

2

-

e

2

n

2

-


(


c

1

n


-

c

2

n



)

2





-
2




e

2

n




(


c

1

n


-

c

2

n



)














(
225
)











R
n



(
0
)


=







z

12

n



cos





α

+


e

3

n



sin






(

α
+

γ
n

+

θ
n


)


+














c

2

n




cos


(

α
+

γ
n


)



+


b

2

n



sin






α
n










(
226
)













where l


sn


is free length of suspension spring, l


s0n


is initial suspension spring deflection at one g, l


wn


is free length of spring component of wheel and l


w0n


is initial wheel spring deflection at one g.




IV. Equations for Entropy Production




Minimum entropy production (for use in the fitness function of the genetic algorithm) is expressed as:












d
β


S

dt

=






-
2





β
.

2

[



α
.



m
b



A
1



A
2


+


m
sn



B
1



{




z
.


6

n




C

αγη





n



-















(


α
.

+


η
.

n


)



z

6

n




S

αγη





n



-


α
.



A
4



}

+


m
an



B
2



{



(


α
.

+


θ
.

n


)



e

1

n




C

αγθ





n



-













α
.



A
6


}

+


m
wn



B
3



{



(


α
.

+


θ
.

n


)



e

3

n




S

αγ





θ





n



-


α
.



A
6



}



]






m
saw2n

+

m
baI

+


m
b



A
1
2


+


m
sn



B
1
2


+


m
an



B
2
2


+


m
wn



B
3
2








(
227
)









d
α


S

dt

=






-
2




α
.

2



{



m
sn



α
.





z
.


6

n




(


z

6

n


+

E

1

n



)



+


m
sn



z

6

n





η
.

n



E

2

n



+











θ
.

n



m
aw1n



H

2

n



}






m
bbI

+

m
sawIn

+


m
sn




z

6

n




(


z

6

n


+

2


E

1

n




)



-

2


m
aw1n



H

1

n









(
228
)









d

η
n



S

dt

=




η
.

n
3


tg






η
n


-


2



η
.

n
2




z
.


6

n





z

6

n


-

d

1

n









(
229
)









d

z

6

n




S

dt

=

2



η
.

n




z
.


6

n

2


tg






η
n






(
230
)









d

z

12

n




S

dt

=



z
.


12

n

2



(


α
.

+


α
.


tg





α

+

2


β
.


tg





β


)






(
231
)













Other Embodiments




Although the foregoing has been a description and illustration of specific embodiments of the invention, various modifications and changes can be made thereto by persons skilled in the art, without departing from the scope and spirit of the invention as defined by the following claims.



Claims
  • 1. A method for controlling a vehicle suspension system, comprising the steps of: measuring first information from said suspension system by using a first plurality of sensors; providing said first information to a first suspension control system, said first suspension control system configured to provide a desired accuracy for said suspension, said first control system providing a first control signal; measuring second information from said suspension system by using a second plurality of sensors, where said second plurality of sensors comprises fewer sensors than said first plurality of sensors; providing said second information to a second suspension control system, said second suspension control system providing a second control signal; and configuring said second suspension control system using said first suspension control signal and said second suspension control signal such that said second suspension control system will provide acceptable control accuracy when used to control said vehicle suspension system without said first suspension control system.
  • 2. The method of claim 1, wherein said step of configuring comprises generating a physical criteria and generating an information criteria.
  • 3. A method for controlling a vehicle suspension system, comprising the steps of: measuring first information from said suspension system by using first plurality of sensors; providing said first information to a first suspension control system, said first suspension control system configured to provide a desired accuracy for said suspension, said first control system providing a first control signal; measuring second information from said suspension system by using a second plurality of sensors, where said second plurality of sensors comprises fewer sensors than said first plurality of sensors; providing said second information to a second suspension control system, said second suspension control system providing a second control signal; and configuring said second suspension control system using said first suspension control signal and said second suspension control signal, wherein said step of configuring comprises generating a physical criteria and generating an information criteria, and wherein said physical criteria is calculated by an entropy model based on thermodynamic properties of said suspension.
  • 4. The method of claim 3, wherein said physical criteria is calculated by an entropy model.
  • 5. The method of claim 4, wherein said second suspension control system is adapted to reduce an entropy production in said suspension.
  • 6. The method of claim 5, wherein said thermodynamic model is based on suspension pitch angle and roll angle.
  • 7. The method of claim 2, wherein said step of configuring further comprises optimizing said second control system using an optimizer, wherein said optimizer uses a genetic algorithm having a fitness function, wherein a portion of said fitness function based on entropy.
  • 8. The method of claim 2, wherein said step of configuring further comprises providing said physical criteria and said information criteria to a genetic algorithm having a fitness function, said fitness function based on entropy.
  • 9. The method of claim 2, wherein said step of configuring further comprises providing a training signal to a fuzzy neural network in said second control system.
  • 10. The method of claim 1, wherein said second plurality of sensors comprises a vertical accelerometer.
  • 11. A method for controlling a vehicle suspension system, comprising the steps of: measuring first information from said suspension system by using first plurality of sensors; providing said first information to a first suspension control system, said first suspension control system configured to provide a desired accuracy for said suspension, said first control system providing a first control signal; measuring second information from said suspension system by using a second plurality of sensors, where said second plurality of sensors comprises fewer sensors than said first plurality of sensors; providing said second information to a second suspension control system, said second suspension control system providing a second control signal; and configuring said second suspension control system using said first suspension control signal and said second suspension control signal, wherein said first plurality of sensors comprises a pitch angle sensor, a roll angle sensor, a position sensor, and an angle sensor.
  • 12. The method of claim 1, wherein said second control signal is provided by a fuzzy neural network, wherein said fuzzy neural network is trained in an off-line mode.
  • 13. A control apparatus configured to control a suspension, said apparatus comprising: suspension control means for generating a suspension control signal based on information from a plurality of sensors measuring said suspension, said suspension control means trained by optimizer means for generating a training signal, said optimizer means generating said training signal using said control signal and an enhanced control signal provided by an enhanced control means, said enhanced control means and said suspension control means used concurrently while training said suspension control means, said suspension control means used without said enhanced control means to control said suspension after said suspension control means has been trained.
  • 14. A control system adapted to control a suspension, comprising: a first plurality of sensors configured to measure first information about said suspension, a first suspension controller configured to receive at least a portion of said first information, said first suspension controller trained in an off-line mode by using said first suspension controller to produce a first control signal and using a second suspension controller to produce a second control signal and developing a training signal from said first control signal and said second control signal, where said first suspension controller is trained to use at least a portion of said first information signal in order to reduce an entropy production of said suspension when said first suspension controller is used in an online mode without said second suspension controller.
  • 15. The apparatus of claim 14, wherein said first suspension controller comprises a fuzzy neural network configured to be trained by a genetic analyzer having a first fitness function, said first fitness function configured to increase mutual information between said first control signal and said second control signal, said second control signal provided by said second controller configured to receive information from a second plurality of sensors, wherein said second plurality of sensors is greater than said first plurality of sensors.
  • 16. The apparatus of claim 15, wherein said genetic analyzer further comprises a second fitness function configured to reduce entropy production rate of said suspension.
  • 17. The apparatus of claim 16, wherein said genetic analyzer is configured to use said second fitness function to realize a node correction in said fuzzy neural network.
  • 18. The apparatus of claim 15, wherein said first plurality of sensors comprises a vertical accelerometer.
  • 19. A control system adapted to control a suspension, comprising: a first plurality of sensors configured to measure first information about said suspension, a first suspension controller configured to receive at least a portion of said first information, said first suspension controller trained to produce a first control signal, where said first suspension controller is trained to use at least a portion of said first information signal in order to reduce an entropy production of said suspension, wherein said first suspension controller comprises a fuzzy neural network configured to be trained by a genetic analyzer having a first fitness function, said first fitness function configured to increase mutual information between said first control signal and a second control signal, said second control signal provided by a second controller configured to receive information from a second plurality of sensors, wherein said second plurality of sensors is greater than said first plurality of sensors, and wherein said second plurality of sensors comprises a pitch angle sensor and a roll angle sensor.
  • 20. A control system adapted to control a suspension, comprising: a first plurality of sensors configured to measure first information about said suspension, a first suspension controller configured to receive at least a portion of said first information, said first suspension controller trained to produce a first control signal, where said first suspension controller is trained to use at least a portion of said first information signal in order to reduce an entropy production of said suspension, wherein said first suspension controller comprises a fuzzy neural network configured to be trained by a genetic analyzer having a first fitness function, said first fitness function configured to increase mutual information between said first control signal and a second control signal, said second control signal provided by a second controller configured to receive information from a second plurality of sensors, wherein said second plurality of sensors is greater than said first plurality of sensors, and wherein said first control signal comprises a hydraulic valve control signal configured to control hydraulic fluid in a shock absorber.
US Referenced Citations (10)
Number Name Date Kind
5208749 Adachi et al. May 1993 A
5324069 Ogawa Jun 1994 A
5483450 Titli Et Al. Jan 1996 A
5488562 Otterbein et al. Jan 1996 A
5557520 Suissa et al. Sep 1996 A
5912821 Kobayashi Jun 1999 A
5928297 Murata et al. Jul 1999 A
5971579 Kim Oct 1999 A
6021369 Kamihira et al. Feb 2000 A
6064996 Yamaguchi et al. May 2000 A
Non-Patent Literature Citations (29)
Entry
Bose Bimal K., (1994), “Expert System, Fuzzy Logic, and Neural Network Applications in Power Electric and Motion Control”, Processing of the IEEE, vol. 82, No. 8, pp. 1303-1323.
Feng Q. and Yamafuji K., (1988), “Design and simulation of control system of an inverted pendulum”, Robotica, vol. 6, No. 3, pp. 235-241.
Gradetsky V.G. and Ulyanov S.V., (1993), “Mobile system with wall climbing robots”, Intern. J. of Computer and Systems Sciences, vol. 31, No. 1, pp. 126-142.
Johansson R., Magnusson M. and Fransson P.A., (1995), “Galvanic vestibular stimulation for analysis of postural adaptation and stability”, IEEE Transactions of Biomedical Engineering, vol. 42, No. 3, pp. 282-292.
Ju M.-S. Yi S.-G,Tsuei Y.-G, and Chou Y.-L., (1995), “Fuzzy control of electrohydraulic above-knee prostheses”, JSME International Journal, vol. 38, No. 1, pp. 78-85.
Lee Y.N., Kim T.W. and Suh I.H., (1994), “A look-up table-based self-organizing fuzzy plus linear controller”, Mechatronics, vol. 4, No. 1, pp. 71-90.
Liu T.S. and Wu J.C., (1993), “A model for rider-motorcycle system using fuzzy control”, IEEE Transactions on Systems, Man and Cybernetics, vol. 23-SMC, No. 1, pp. 267-276.
Mendel, Jerry M., (1995), “Fuzzy Logic Systems for Engineering: A Tutorial”, Proceeding of the IEEE, vol. 83, No. 3, pp. 345-377.
Nakajima R., Tsubouchi T., Yuta S. and Koyanagi E., (1997), “A development of a new mechanism of an autonomous unicycle”, Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS '97), vol. 2, Grenoble, France, pp. 906-912.
Perround M. and Saucier A., (1987), “Thermodynamics of dissipative systems”, Helvetica Physica, vol. 60, No. 8, pp. 1038-1051.
Sheng Z.Q. and Yamafuji K., (1995), “Realization of a human riding a unicycle by a robot”, Proc. of '95 IEEE Intern. Conf. on Robotics & Automation, Japan, vol. 2, pp. 1319-1326.
Sheng Z.Q. and Yamafuji K., (1997), “Postural stability of a human riding a unicycle and its emulation by a robot”, IEEE Trans. on Robotics and Automation, vol. 13, No. 5, pp. 709-720.
Sheng, Z.Q. and Yamafuji K., (1995), “Study on the stability and motion control of a unicycle, 1st Report: Dynamics of a human riding a unicycle and its modeling by link mechanism”, JSME International Journal, vol. 38C, No. 2, pp. 249-259.
Sheng, Z.Q., Yamafuji K. and Ulyanov S.V., (1996), “Study on the stability and motion control of a unicycle. Pts 3,” JSME International Journal, vol. 39, No. 3, pp. 560-568.
Tsai, C.H. and Li, T.H.S., (1996) “Design and Lyapunov Function Based Fuzzy Logic Controller”, IEEE, pp.396-401.
Ulyanov S.V., (1992), “Fuzzy models of Intelligent control systems: Theoretical and applied aspects” Soviet Journal of Computer and Systems Sciences (Engineering Cybernetics), vol. 30, No. 4, pp. 1-22.
Ulyanov S.V., Sheng Z.Q., Yamafuji K., Watanabe S. and Ohkura T., (1995), “Self-organization fuzzy chaos Intelligent controller for a robotic unicycle: A New benchmark in AI control”, Proc. of 5th Intelligent System Symposium: Fuzzy, AI and Neural Network Applications Technologies (FAN Symp. '95), Tokyo, pp. 41-46.
Ulyanov S.V. and Yamafuji K., (1996), “Fuzzy Intelligent emotion and instinct control of a robotic unicyle”, 4th Intern. Workshop on Advanced Motion Control, Mie, Japan, vol. 1, pp. 127-132.
Ulyanov S.V., Feng Q., Yamafuji K. and Ulyanov V.S., (1998), “Stochastic analysis of nonlinear dynamic system with time-dependent structure, Pts. 1,2”, Probabilistic Engineering Mechanics (published).
Ulyanov, S.V., Yamafuji, K., Gradetsky V.G., and Pagni, A., “Expert Fuzzy-Neuro Controller Desing for Wall Climbing Robot for Decontamination of Nuclear-Power Station”, Journal of Robotics and Mechatronics, vol. 7, No. 1, 1995, pp. 75-86.
V.S. Ulyanov, Watanabe, S. Yamafuji, K., Ulyanov, S.V., Litvintseva, L.V., and Kurawaki, I. “Intelligent robust control of a robotic unicycle based on a new physical measure for mechanical controllability” Advanced Robotics, vol. 12, No. 4 1998, pp. 455-481.
Ulyanov, S.V., Yamafuji, K., Miyagawa, K., Tanaka, T. and Fukuda T., “Intelligent Fuzzy Motion Control of Mobile Robot for Service Use”, IEEE/RSJ International Conference on Intelligent Robots and Systems, Human Robot Interaction and Cooperative Robots, vol. 3, Aug. 5-9, 1995 pp.486-491.
Ulyanov, S.V., Yamafuji, K., Arai, F., and Fukuda T., “Modelling of Micro-Nano-Robots and Physical Limit of Micro Control” Journal of the Robotics Society of Japan, vol. 14, No. 8, Nov. 1996, pp. 18-21.
van Rooij A.J.F., Jain L.C. and Johnson R.P., (1997), “Neural Network Training Using Genetic Algorithms” (Machine Perception and Artificial Intelligence, vol. 26). World Scientific Pub., SGP.
Vasilieva O.I., Ionov I.P. and Ulyanov S.V., (1989), “Dual control of the artificial ventilation process with use of a fuzzy controller in the feedback circuit”, Biomedical Eng., vol. 23, No. 1, pp. 7-16.
Wu J.C. and Llu T.S., (1995), “Fuzzy control of rider-motorcycle system using genetic algorithm and auto-tuning”, Mechatronics, vol. 4, No. 4, pp. 441-455.
Wu J.C. and Llu T.S., (1996), “Fuzzy control stabilization with applications to motorcycle control”, IEEE Transactions on Systems, Man and Cybernetics, vol. 26, No. 6, pp. 836-847.
Zhakharov V.N. and Ulyanov S.V., (1994), “Fuzzy models of intelligent industrial controllers and controls systems. Pts 1”, Journal of Computer and Systems Sciences International, vol. 32, No. 1, pp. 123-144.
Zhakharov V.N. and Ulyanov S.V., (1995), “Fuzzy models of intelligent industrial controllers and control systems. Pts 2,3”, Journal of Computer and Systems Sciences International, vol. 33, No. 2, pp. 94-108; 117-136.