Patent Publication Number: US-6668214-B2

Title: Control system for throttle valve actuating device

Description:
RELATED APPLICATIONS 
     This application is a divisional of U.S. patent application Ser. No. 10/043,625, filed on Jan. 10, 2002, and entitled “CONTROL SYSTEM FOR PLANT” now abandoned. 
    
    
     BACKGROUND OF THE INVENTION 
     The present invention relates to a control system for a throttle valve actuating device including a throttle valve of an internal combustion engine and an actuator for actuating the throttle valve. 
     One known control system for controlling a throttle valve actuating device including a motor for actuating the throttle valve and springs for energizing the throttle valve is disclosed in Japanese Patent Laid-open No. Hei 8-261050. In this control system, the throttle valve actuating device is controlled with a PID (Proportional, Integral, and Differential) control according to a throttle valve opening detected by a throttle valve opening sensor. 
     Since the throttle valve actuating device includes a dead time element which delays the throttle valve movement, the PID control by the above-described conventional control system may easily become unstable if the control gain of the PID control according to the detected throttle valve opening is set to a large value. Therefore, it is necessary to set the control gain to a value which is sufficiently small so that the control may not become unstable. However, this results in a low response speed. 
     SUMMARY OF THE INVENTION 
     It is therefore an object of the present invention to provide a control system for a throttle valve actuating device, which can improve the response speed by expanding the stability margin of the control. 
     To achieve the above object, the present invention provides a control system for a throttle valve actuating device ( 10 ) including a throttle valve ( 3 ) of an internal combustion engine and actuating means ( 6 ) for actuating the throttle valve ( 3 ). The control system includes predicting means ( 23 ) for predicting a future throttle valve opening (PREDTH) and controls the throttle actuating device ( 10 ) according to the throttle valve opening (PREDTH) predicted by the predicting means ( 23 ) so that the throttle valve opening (TH) coincides with a target opening (THR). 
     With this configuration, the throttle actuating device is controlled according to the future throttle valve opening predicted by the predicting means. Accordingly, the stability margin of the control can be expanded although the throttle valve actuating device includes a dead time element. As a result, the control gain can be set to a larger value to thereby improve the control response speed. 
     Preferably, the throttle valve actuating device ( 10 ) is modeled to a controlled object model which includes a dead time element, and the predicting means ( 23 ) predicts a throttle valve opening after the elapse of a dead time (d) due to the dead time element, based on the controlled object model. 
     With this configuration, the throttle valve actuating device ( 10 ) is modeled to the controlled object model including a dead time element, and the throttle valve opening after the elapse of a dead time is predicted based on the controlled object model. Accordingly, compensation of the dead time in the throttle valve actuating device can accurately be performed. 
     Preferably, the control system further includes identifying means ( 22 ) for identifying at least one model parameter (a 1 , a 2 , b 1 , c 1 ) of the controlled object model. The predicting means ( 23 ) predicts the throttle valve opening (PREDTH) using the at least one model parameter (a 1 , a 2 , b 1 , c 1 ) identified by the identifying means ( 22 ). 
     With this configuration, the throttle valve opening is predicted using one or more model parameter identified by the identifying means. Accordingly, it is possible to calculate an accurate predicted throttle valve opening even when the dynamic characteristic of the throttle valve actuating device has changed due to aging or a change in the environmental condition. 
     Preferably, the control system further includes a sliding mode controller ( 21 ) for controlling the throttle valve actuating device ( 10 ) with a sliding mode control according to a throttle valve opening (PREDTH) predicted by the predicting means ( 23 ). 
     With this configuration, the throttle valve actuating device is controlled with the sliding mode control having good robustness. Accordingly, good stability and controllability of the control can be maintained even when there exists a modeling error (a difference between the characteristics of the actual plant and the characteristics of the controlled object model) due to a difference between the actual dead time of the throttle valve actuating device and the dead time of the controlled object model. 
     Preferably, the sliding mode controller ( 21 ) controls the throttle valve actuating device ( 10 ) using the at least one model parameter (a 1 , a 2 , b 1 , c 1 ) identified by the identifying means ( 22 ). 
     With this configuration, the throttle valve actuating device is controlled by the sliding mode controller using one or more model parameter identified by the identifying means. Accordingly, the modeling error can be made smaller so that the controllability is further improved. 
     Preferably, the control input (Usl) from the sliding mode controller ( 21 ) to the throttle valve actuating device ( 10 ) includes an adaptive law input (Uadp). 
     With this configuration, better controllability is obtained even in the presence of disturbance and/or the modeling error. 
     The above and other objects, features, and advantages of the present invention will become apparent from the following description when taken in conjunction with the accompanying drawings which illustrate embodiments of the present invention by way of example. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a schematic view of a throttle valve control system according to a first embodiment of the present invention; 
     FIGS. 2A and 2B are diagrams showing frequency characteristics of the throttle valve actuating device shown in FIG. 1; 
     FIG. 3 is a functional block diagram showing functions realized by an electronic control unit (ECU) shown in FIG. 1; 
     FIG. 4 is a diagram showing the relationship between control characteristics of a sliding mode controller and the value of a switching function setting parameter (VPOLE); 
     FIG. 5 is a diagram showing a range for setting control gains (F, G) of the sliding mode controller; 
     FIGS. 6A and 6B are diagrams illustrative of a drift of model parameters; 
     FIGS. 7A through 7C are diagrams showing functions for correcting an identifying error; 
     FIG. 8 is a diagram illustrating that a default opening deviation of a throttle valve is reflected to a model parameter (c 1 ′); 
     FIG. 9 is a flowchart showing a throttle valve opening control process; 
     FIG. 10 is a flowchart showing a process of setting state variables in the process shown in FIG. 9; 
     FIG. 11 is a flowchart showing a process of performing calculations of a model parameter identifier in the process shown in FIG. 9; 
     FIG. 12 is a flowchart showing a process of calculating an identifying error (ide) in the process shown in FIG. 11; 
     FIGS. 13A and 13B are diagrams illustrative of a process of low-pass filtering on the identifying error (ide); 
     FIG. 14 is a flowchart showing the dead zone process in the process shown in FIG. 12; 
     FIG. 15 is a diagram showing a table used in the process shown in FIG. 14; 
     FIG. 16 is a flowchart showing a process of stabilizing a model parameter vector (θ) in the process shown in FIG. 11; 
     FIG. 17 is a flowchart showing a limit process of model parameters (a 1 ′, a 2 ′) in the process shown in FIG. 16; 
     FIG. 18 is a diagram illustrative of the change in the values of the model parameters in the process shown in FIG. 16; 
     FIG. 19 is a flowchart showing a limit process of a model parameter (b 1 ′) in the process shown in FIG. 16; 
     FIG. 20 is a flowchart showing a limit process of a model parameter (c 1 ′) in the process shown in FIG. 16; 
     FIG. 21 is a flowchart showing a process of performing calculations of a state predictor in the process shown in FIG. 9; 
     FIG. 22 is a flowchart showing a process of calculating a control input (Usl) in the process shown in FIG. 9; 
     FIG. 23 is a flowchart showing a process of calculating a predicted switching function value (σpre) in the process shown in FIG. 22; 
     FIG. 24 is a flowchart showing a process of calculating the switching function setting parameter (VPOLE) in the process shown in FIG. 23; 
     FIGS. 25A through 25C are diagrams showing maps used in the process shown in FIG. 24; 
     FIG. 26 is a flowchart showing a process of calculating an integrated value of the predicted switching function value (σpre) in the process shown in FIG. 22; 
     FIG. 27 is a flowchart showing a process of calculating a reaching law input (Urch) in the process shown in FIG. 22; 
     FIG. 28 is a flowchart showing a process of calculating an adaptive law input (Uadp) in the process shown in FIG. 22; 
     FIG. 29 is a flowchart showing a process of determining the stability of the sliding mode controller in the process shown in FIG. 9; 
     FIG. 30 is a flowchart showing a process of calculating a default opening deviation (thdefadp) in the process shown in FIG.  9 ; 
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     FIG. 1 schematically shows a configuration of a throttle valve control system according to a first embodiment of the present invention. An internal combustion engine (hereinafter referred to as “engine”)  1  has an intake passage  2  with a throttle valve  3  disposed therein. The throttle valve  3  is provided with a return spring  4  as a first energizing means for energizing the throttle valve  3  in a closing direction, and a resilient member  5  as a second energizing means for energizing the throttle valve  3  in an opening direction. The throttle valve  3  can be actuated by a motor  6  as an actuating means through gears (not shown). When the actuating force from the motor  6  is not applied to the throttle valve  3 , an opening TH of the throttle valve  3  is maintained at a default opening THDEF (for example, 5 degrees) where the energizing force of the return spring  4  and the energizing force of the resilient member  5  are in equilibrium. 
     The motor  6  is connected to an electronic control unit (hereinafter referred to as “ECU”)  7 . The operation of the motor  6  is controlled by the ECU  7 . The throttle valve  3  is associated with a throttle valve opening sensor  8  for detecting the throttle valve opening TH. A detected signal from the throttle valve opening sensor  8  is supplied to the ECU  7 . 
     Further, the ECU  7  is connected to an acceleration sensor  9  for detecting a depression amount ACC of an accelerator pedal to detect an output demanded by the driver of the vehicle on which the engine  1  is mounted. A detected signal from the acceleration sensor  9  is supplied to the ECU  7 . 
     The ECU  7  has an input circuit, an A/D converter, a central processing unit (CPU), a memory circuit, and an output circuit. The input circuit is supplied with detected signals from the throttle valve opening sensor  8  and the acceleration sensor  9 . The A/D converter converts input signals into digital signals. The CPU carries out various process operations. The memory circuit has a ROM (read only memory) for storing processes executed by the CPU, and maps and tables that are referred to in the processes, a RAM for storing results of executing processes by the CPU. The output circuit supplies an energizing current to the motor  6 . The ECU  7  determines a target opening THR of the throttle valve  3  according to the depression amount ACC of the accelerator pedal, determines a control quantity DUT for the motor  6  in order to make the detected throttle valve opening TH coincide with the target opening THR, and supplies an electric signal according to the control quantity DUT to the motor  6 . 
     In the present embodiment, a throttle valve actuating device  10  that includes the throttle valve  3 , the return spring  4 , the resilient member  5 , and the motor  6  is a controlled object. An input to be applied to the controlled object is a duty ratio DUT of the electric signal applied to the motor  6 . An output from the controlled object is the throttle valve opening TH detected by the throttle valve opening sensor  8 . 
     When frequency response characteristics of the throttle valve actuating device  10  are measured, gain characteristics and phase characteristics indicated by the solid lines in FIGS. 2A and 2B are obtained. A model defined by the equation (1) shown below is set as a controlled object model. Frequency response characteristics of the model are indicated by the broken-line curves in FIGS. 2A and 2B. It has been confirmed that the frequency response characteristics of the model are similar to the characteristics of the throttle valve actuating device  10 . 
     
       
           DTH ( k+ 1)= a   1 × DTH ( k )+ a   2 × DTH ( k− 1)+ b   1 × DUT ( k−d )+ c   1   (1)  
       
     
     where k is a parameter representing discrete time, and DTH(k) is a throttle valve opening deviation amount defined by the equation (2) shown below. DTH(k+1) is a throttle valve opening deviation amount at a discrete time (k+1). 
     
       
           DTH ( k )= TH ( k )− THDEF   (2)  
       
     
     where TH is a detected throttle valve opening, and THDEF is the default opening. 
     In the equation (1), a 1 , a 2 , b 1 , and c 1  are parameters determining the characteristics of the controlled object model, and d is a dead time. The dead time is a delay between the input and output of the controlled object model. 
     The model defined by the equation (1) is a DARX model (delayed autoregressive model with exogeneous input) of a discrete time system, which is employed for facilitating the application of an adaptive control. 
     In the equation (1), the model parameter c 1  which is irrelevant to the input and output of the controlled object is employed, in addition to the model parameters a 1  and a 2  which are relevant to the output deviation amount DTH and the model parameter b 1  which is relevant to the input duty ratio DUT. The model parameter c 1  is a parameter representing a deviation amount of the default opening THDEF and disturbance applied to the throttle valve actuating device  10 . In other words, the default opening deviation amount and the disturbance can be identified by identifying the model parameter c 1  simultaneously with the model parameters a 1 , a 2 , and b 1  by a model parameter identifier. 
     FIG. 3 is a functional block diagram of the throttle valve control system which is realized by the ECU  7 . The throttle valve control system as configured includes, an adaptive sliding mode controller  21 , a model parameter identifier  22 , a state predictor  23  for calculating a predicted throttle valve opening deviation amount (hereinafter referred to as “predicted deviation amount” or PREDTH(k)) where PREDTH(k) (=DTH(k+d)) after the dead time d has elapsed, and a target opening setting unit  24  for setting a target opening THR for the throttle valve  3  according to the accelerator pedal depression amount ACC. 
     The adaptive sliding mode controller  21  calculates a duty ratio DUT according to an adaptive sliding mode control in order to make the detected throttle valve opening TH coincide with the target opening THR, and outputs the calculated duty ratio DUT. 
     By using the adaptive sliding mode controller  21 , it is possible to change the response characteristics of the throttle valve opening TH to the target opening THR, using a specific parameter (VPOLE). As a result, it is possible to avoid shocks at the time the throttle valve  3  moves from an open position to a fully closed position, i.e., at the time the throttle valve  3  collides with a stopper for stopping the throttle valve  3  in the fully closed position. It is also possible to make the engine response corresponding to the operation of the accelerator pedal variable. Further, it is also possible to obtain a good stability against errors of the model parameters. 
     The model parameter identifier  22  calculates a corrected model parameter vector θL (θL T =[a 1 , a 2 , b 1 , c 1 ]) and supplies the calculated corrected model parameter vector θL to the adaptive sliding mode controller  21 . More specifically, the model parameter identifier  22  calculates a model parameter vector θ based on the throttle valve opening TH and the duty ratio DUT. The model parameter identifier  22  then carries out a limit process of the model parameter vector θ to calculate the corrected model parameter vector θL, and supplies the corrected model parameter vector θL to the adaptive sliding mode controller  21 . In this manner, the model parameters a 1 , a 2 , and b 1  which are optimum for making the throttle valve opening TH follow up the target opening THR are obtained, and also the model parameter c 1  indicative of disturbance and a deviation amount of the default opening THDEF is obtained. 
     By using the model parameter identifier  22  for identifying the model parameters on a real-time basis, adaptation to changes in engine operating conditions, compensation for hardware characteristics variations, compensatation for power supply voltage fluctuations, and adaptation to aging-dependent changes of hardware characteristics are possible. 
     The state predictor  23  calculates a throttle valve opening TH (predicted value) after the dead time d has elapsed, or more specifically a predicted deviation amount PREDTH, based on the throttle valve opening TH and the duty ratio DUT, and supplies the calculated deviation amount PREDTH to the adaptive sliding mode controller  21 . By using the predicted deviation amount PREDTH, the robustness of the control system against the dead time of the controlled object is ensured, and the controllability in the vicinity of the default opening THDEF where the dead time is large is improved. 
     Next, principles of operation of the adaptive sliding mode controller  21  will be hereinafter described. 
     First, a target value DTHR(k) is defined as a deviation amount between the target opening THR(k) and the default opening THDEF by the following equation (3). 
     
       
           DTHR ( k )= THR ( k )− THDEF   (3)  
       
     
     If a deviation e(k) between the throttle valve opening deviation amount DTH and the target value DTHR is defined by the following equation (4), then a switching function value σ(k) of the adaptive sliding mode controller is set by the following equation (5).                e        (   k   )       =       DTH        (   k   )       -     DTHR        (   k   )                 (   4   )                       σ        (   k   )       =                  e        (   k   )       +     VPOLE   ×     e        (     k   -   1     )                       =                  DTH        (   k   )       -     DTHR        (   k   )       +                              VPOLE   ×     (       DTH        (     k   -   1     )       -     DTHR        (     k   -   1     )         )                     (   5   )                         
     where VPOLE is a switching function setting parameter that is set to a value which is greater than −1 and less than 1. 
     On a phase plane defined by a vertical axis representing the deviation e(k) and a horizontal axis representing the preceding deviation e(k−1), a pair of the deviation e(k) and the preceding deviation e(k−1) satisfying the equation of “σ(k)=0” represents a straight line. The straight line is generally referred to as a switching straight line. A sliding mode control is a control contemplating the behavior of the deviation e(k) on the switching straight line. The sliding mode control is carried out so that the switching function value σ(k) becomes 0, i.e., the pair of the deviation e(k) and the preceding deviation e(k−1) exists on the switching straight line on the phase plane, to thereby achieve a robust control against disturbance and the modeling error (the difference between the characteristics of an actual plant and the characteristics of a controlled object model). As a result, the throttle valve opening deviation amount DTH is controlled with good robustness to follow up the target value DTHR. 
     As shown in FIG. 4, by changing the value of the switching function setting parameter VPOLE in the equation (5), it is possible to change damping characteristics of the deviation e(k), i.e., the follow-up characteristics of the throttle valve opening deviation amount DTH to follow the target value DTHR. Specifically, if VPOLE equals −1, then the throttle valve opening deviation amount DTH completely fails to follow up the target value DTHR. As the absolute value of the switching function setting parameter VPOLE is reduced, the speed at which the throttle valve opening deviation amount DTH follows up the target value DTHR increases. 
     The throttle valve control system is required to satisfy the following requirements A1 and A2: 
     A1) When the throttle valve  3  is shifted to the fully closed position, collision of the throttle valve  3  with the stopper for stopping the throttle valve  3  in the fully closed position should be avoided; and 
     A2) The controllability with respect to the nonlinear characteristics in the vicinity of the default opening THDEF (a change in the resiliency characteristics due to the equilibrium between the energizing force of the return spring  4  and the energizing force of the resilient member  5 , backlash of gears interposed between the motor  6  and the throttle valve  3 , and a dead zone where the throttle valve opening does not change even when the duty ratio DUT changes) should be improved. 
     Therefore, it is necessary to lower the speed at which the deviation e(k) converges, i.e., the converging speed of the deviation e(k), in the vicinity of the fully closed position of the throttle valve, and to increase the converging speed of the deviation e(k) in the vicinity of the default opening THDEF. 
     According to the sliding mode control, the converging speed of e(k) can easily be changed by changing the switching function setting parameter VPOLE. Therefore in the present embodiment, the switching function setting parameter VPOLE is set according to the throttle valve opening TH and an amount of change DDTHR (=DTHR(k)−DTHR(k−1)) of the target value DTHR, to thereby satisfy the requirements A1 and A2. 
     As described above, according to the sliding mode control, the deviation e(k) is converged to 0 at an indicated converging speed and robustly against disturbance and the modeling error by constraining the pair of the deviation e(k) and the preceding deviation e(k−1) on the switching straight line (the pair of e(k) and e(k−1) will be hereinafter referred to as “deviation state quantity”). Therefore, in the sliding mode control, it is important how to place the deviation state quantity onto the switching straight line and constrain the deviation state quantity on the switching straight line. 
     From the above standpoint, an input DUT(k) (also indicated as Usl(k)) to the controlled object (an output of the controller) is expressed as the sum of an equivalent control input Ueq(k), a reaching law input Urch(k), and an adaptive law input Uadp(k), as indicated by the following equation (6).                      DUT        (   k   )       =                Usl        (   k   )                   =                  Ueq        (   k   )       +     Urch        (   k   )       +     Uadp        (   k   )                       (   6   )                         
     The equivalent control input Ueq(k) is an input for constraining the deviation state quantity on the switching straight line. The reaching law input Urch(k) is an input for placing deviation state quantity onto the switching straight line. The adaptive law input Uadp(k) is an input for placing deviation state quantity onto the switching straight line while reducing the modeling error and the effect of disturbance. Methods of calculating these inputs Ueq(k), Urch(k), and Uadp(k) will be described below. 
     Since the equivalent control input Ueq(k) is an input for constraining the deviation state quantity on the switching straight line, a condition to be satisfied is given by the following equation (7): 
     
       
         σ( k )=σ( k+ 1)  (7)  
       
     
     Using the equations (1), (4), and (5), the duty ratio DUT(k) satisfying the equation (7) is determined by the equation (9) shown below. The duty ratio DUT(k) calculated with the equation (9) represents the equivalent control input Ueq(k). The reaching law input Urch(k) and the adaptive law input Uadp(k) are defined by the respective equations (10) and (11) shown below.                      DUT        (   k   )       =                  1   b1          {         [     1   -   a1   -   VPOLE     )          DTH        (     k   +   d     )         +                                      (     VPOLE   -   a2     )          DTH        (     k   +   d   -   1     )         -   c1   +                                DTHR        (     k   +   d   +   1     )       +       (     VPOLE   -   1     )          DTHR        (     k   +   d     )         -                              VPOLE   ×     DTHR        (     k   +   d   -   1     )         }               =                Ueq        (   k   )                     (   9   )                 Urch        (   k   )       =         -   F     b1                     σ        (     k   +   d     )                 (   10   )                 Uadp        (   k   )       =         -   G     b1            ∑     i   =   0       k   +   d            Δ                 T                   σ        (   i   )                     (   11   )                         
     where F and G respectively represent a reaching law control gain and an adaptive law control gain, which are set as described below, and ΔT represents a control period. 
     Calculating the equation (9) requires a throttle valve opening deviation amount DTH(k+d) after the elapse of the dead time d and a corresponding target value DTHR(k+d+1). Therefore, the predicted deviation amount PREDTH(k) calculated by the state predictor  23  is used as the throttle valve opening deviation amount DTH(k+d) after the elapse of the dead time d, and the latest target value DTHR is used as the target value DTHR(k+d+1). 
     Next, the reaching law control gain F and the adaptive law control gain G are determined so that the deviation state quantity can stably be placed onto the switching straight line by the reaching law input Urch and the adaptive law input Uadp. 
     Specifically, a disturbance V(k) is assumed, and a stability condition for keeping the switching function value σ(k) stable against the disturbance V(k) are determined to obtain a condition for setting the gains F and G. As a result, it has been obtained as the stability condition that the combination of the gains F and G satisfies the following equations (12) through (14), in other words, the combination of the gains F and G should be located in a hatched region shown in FIG.  5 . 
     
       
         F&gt;0  (12)  
       
     
     
       
         G&gt;0  (13)  
       
     
     
       
           F&lt; 2−(Δ T/ 2) G   (14)  
       
     
     As described above, the equivalent control input Ueq(k), the reaching law input Urch(k), and the adaptive law input Uadp(k) are calculated from the equations (9) through (11), and the duty ratio DUT(k) is calculated as the sum of those inputs. 
     The model parameter identifier  22  calculates a model parameter vector of the controlled object model, based on the input (DUT(k)) and output (TH(k)) of the controlled object, as described above. Specifically, the model parameter identifier  22  calculates a model parameter vector θ(k) according to a sequential identifying algorithm (generalized sequential method-of-least-squares algorithm) represented by the following equation (15). 
     
       
         θ( k )=θ( k− 1)+ KP ( k ) ide ( k )  (15)  
       
     
     
       
         θ( k ) T   =[a   1 ′,  a   2 ′,  b   1 ′,  c   1 ′]  (16)  
       
     
     where a 1 ′, a 2 ′, b 1 ′, c 1 ′ represent model parameters before a limit process described later is carried out, ide(k) represents an identifying error defined by the equations (17), (18), and (19) shown below, where DTHHAT(k) represents an estimated value of the throttle valve opening deviation amount DTH(k) (hereinafter referred to as “estimated throttle valve opening deviation amount”) which is calculated using the latest model parameter vector θ(k−1), and KP(k) represents a gain coefficient vector defined by the equation (20) shown below. In the equation (20), P(k) represents a quartic square matrix calculated from the equation (21) shown below. 
     
       
           i d e  ( k )= DTH ( k )− DTHHAT ( k )  (17)  
       
     
     
       
           DTHHAT ( k )=θ( k− 1) T ζ( k )  (18)  
       
     
     
       
         ζ( k ) T   =[DTH ( k− 1),  DTH ( k− 2),  DUT ( k−d− 1), 1]  (19)  
       
       
         
           
             
               
                 
                   
                     K 
                      
                     
                         
                     
                      
                     
                       P 
                        
                       
                         ( 
                         k 
                         ) 
                       
                     
                   
                   = 
                   
                     
                       
                         P 
                          
                         
                           ( 
                           k 
                           ) 
                         
                       
                        
                       
                         ζ 
                          
                         
                           ( 
                           k 
                           ) 
                         
                       
                     
                     
                       1 
                       + 
                       
                         
                           
                             ζ 
                             T 
                           
                            
                           
                             ( 
                             k 
                             ) 
                           
                         
                          
                         
                           P 
                            
                           
                             ( 
                             k 
                             ) 
                           
                         
                          
                         
                           ζ 
                            
                           
                             ( 
                             k 
                             ) 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   20 
                   ) 
                 
               
             
             
               
                 
                   
                     P 
                      
                     
                       ( 
                       
                         k 
                         + 
                         1 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       1 
                       
                         λ 
                         1 
                       
                     
                      
                     
                       ( 
                       
                         - 
                         
                           
                             
                               λ 
                               2 
                             
                              
                             
                               P 
                                
                               
                                 ( 
                                 k 
                                 ) 
                               
                             
                              
                             
                               ζ 
                                
                               
                                 ( 
                                 k 
                                 ) 
                               
                             
                              
                             
                               
                                 ζ 
                                 T 
                               
                                
                               
                                 ( 
                                 k 
                                 ) 
                               
                             
                           
                           
                             
                               λ 
                               1 
                             
                             + 
                             
                               
                                 λ 
                                 2 
                               
                                
                               
                                 
                                   ζ 
                                   T 
                                 
                                  
                                 
                                   ( 
                                   k 
                                   ) 
                                 
                               
                                
                               
                                 P 
                                  
                                 
                                   ( 
                                   k 
                                   ) 
                                 
                               
                                
                               
                                 ζ 
                                  
                                 
                                   ( 
                                   k 
                                   ) 
                                 
                               
                             
                           
                         
                       
                       ) 
                     
                      
                     
                       P 
                        
                       
                         ( 
                         k 
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
             
               
                 
                   ( 
                   
                     : 
                     
                       I 
                        
                       
                           
                       
                        
                       d 
                        
                       
                           
                       
                        
                       e 
                        
                       
                           
                       
                        
                       n 
                        
                       
                           
                       
                        
                       t 
                        
                       
                           
                       
                        
                       i 
                        
                       
                           
                       
                        
                       t 
                        
                       
                           
                       
                        
                       y 
                        
                       
                           
                       
                        
                       M 
                        
                       
                           
                       
                        
                       a 
                        
                       
                           
                       
                        
                       t 
                        
                       
                           
                       
                        
                       r 
                        
                       
                           
                       
                        
                       i 
                        
                       
                           
                       
                        
                       x 
                     
                   
                   ) 
                 
               
               
                 
                     
                 
               
             
           
         
         
         
             
         
       
     
     In accordance with the setting of coefficients λ1 and λ2 in the equation (21), the identifying algorithm from the equations (15) through (21) becomes one of the following four identifying algorithm: 
     
       
         
           
               
               
             
               
                   
               
             
            
               
                 λ1 = 1, λ2 = 0 
                 Fixed gain algorithm 
               
               
                 λ1 = 1, λ2 = 1 
                 Method-of-least-squares 
               
               
                   
                 algorithm 
               
               
                 λ1 = 1, λ2 = λ 
                 Degressive gain 
               
               
                   
                 algorithm 
               
               
                   
                 (λ is a given value 
               
               
                   
                 other than 0, 1) 
               
               
                 λ1 = λ, λ2 = 1 
                 Weighted Method-of- 
               
               
                   
                 least- squares algorithm 
               
               
                   
                 (λ is a given value 
               
               
                   
                 other than 0, 1) 
               
               
                   
               
            
           
         
       
     
     In the present embodiment, it is required that the following requirements B1, B2, and B3 are satisfied: 
     B1) Adaptation to quasi-static dynamic characteristics changes and hardware characteristics variations 
     “Quasi-static dynamic characteristics changes” mean slow-rate characteristics changes such as power supply voltage fluctuations or hardware degradations due to aging. 
     B2) Adaptation to high-rate dynamic characteristics changes 
     Specifically, this means adaptation to dynamic characteristics changes depending on changes in the throttle valve opening TH. 
     B3) Prevention of a drift of model parameters 
     The drift, which is an excessive increase of the absolute values of model parameters, should be prevented. The drift of model parameters is caused by the effect of the identifying error, which should not be reflected to the model parameters, due to nonlinear characteristics of the controlled object. 
     In order to satisfy the requirements B1 and B2, the coefficients λ1 and λ2 are set respectively to a given value λ and “0” so that the weighted Method-of-least-squares algorithm is employed. 
     Next, the drift of model parameters will be described below. As shown in FIG.  6 A and FIG. 6B, if a residual identifying error, which is caused by nonlinear characteristics such as friction characteristics of the throttle valve, exists after the model parameters have been converged to a certain extent, or if a disturbance whose average value is not zero is steadily applied, then residual identifying errors are accumulated, causing a drift of model parameters. 
     Since such a residual identifying error should not be reflected to the values of model parameters, a dead zone process is carried out using a dead zone function Fn 1  as shown in FIG.  7 A. Specifically, a corrected identifying error idenl(k) is calculated from the following equation (23), and a model parameter vector θ(k) is calculated using the corrected identifying error idenl(k). That is, the following equation (15a) is used instead of the above equation (15). In this manner, the requirement B3) can be satisfied. 
     
       
           idenl ( k )= Fnl ( ide ( k ))  (23)  
       
     
     
       
         θ( k )=θ( k− 1)+ KP ( k ) idenl ( k )  (15a)  
       
     
     The dead zone function Fn 1  is not limited to the function shown in FIG. 7A. A discontinuous dead zone function as shown in FIG. 7B or an incomplete dead zone function as shown in FIG. 7C may be used as the dead zone function Fnl. However, it is impossible to completely prevent the drift if the incomplete dead zone function is used. 
     The amplitude of the residual identifying error changes according to an amount of change in the throttle valve opening TH. In the present embodiment, a dead zone width parameter EIDNRLMT which defines the width of the dead zone shown in FIGS. 7A through 7C is set according to the square average value DDTHRSQA of an amount of change in the target throttle valve opening THR. Specifically, the dead zone width parameter EIDNRLMT is set such that it increases as the square average value DDTHRSQA increases. According to such setting of the dead zone width parameter EIDNRLMT, it is prevented to neglect an identifying error to be reflected to the values of the model parameters as the residual identifying error. In the following equation (24), DDTHR represents an amount of change in the target throttle valve opening THR, which is calculated from the following equation (25):                D                 D                 T                 H                 R                 S                 Q                   A        (   k   )         =       1     n   +   1              ∑     i   =   0     n          D                 D                 T                 H                     R        (   i   )       2                   (   24   )                         D                 D                 T                 H                   R        (   k   )         =       D                 T                 H                   R        (   k   )         -     D                 T                 H                   R        (     k   -   1     )                       =       T                 H                   R        (   k   )         -     T                 H                   R        (     k   -   1     )                              (   25   )                         
     Since the throttle valve opening deviation amount DTH is controlled to the target value DTHR by the adaptive sliding mode controller  21 , the target value DTHR in the equation (25) may be changed to the throttle valve opening deviation amount DTH. In this case, an amount of change DDTH in the throttle valve opening deviation amount DTH may be calculated, and the dead zone width parameter EIDNRLMT may be set according to the square average value DDTHRSQA obtained by replacing DDTHR in the equation (24) with DDTH. 
     For further improving the robustness of the control system, it is effective to further stabilize the adaptive sliding mode controller  21 . In the present embodiment, the elements a 1 ′, a 2 ′, b 1 ′, and c 1 ′ of the model parameter vector θ(k) calculated from the equation (15) are subjected to the limit process so that a corrected model parameter vector θL(k)(θL(k) T =[a 1 , a 2 , b 1 , c 1 ]) is calculated. The adaptive sliding mode controller  21  performs a sliding mode control using the corrected model parameter vector θL(k). The limit process will be described later in detail referring to the flowcharts. 
     Next, a method for calculating the predicted deviation amount PREDTH in the state predictor  23  will be described below. 
     First, matrixes A, B and vectors X(k), U(k) are defined according to the following equations (26) through (29).              A   =     [         a1       a2           1       0         ]             (   26   )               B   =     [         b1       c1           0       0         ]             (   27   )                 X        (   k   )       =     [           D                 T                   H        (   k   )                   D                 T                   H        (     k   -   1     )               ]             (   28   )                 U        (   k   )       =     [           D                 U                   T        (   k   )                 1         ]             (   29   )                         
     By rewriting the equation (1) which defines the controlled object model, using the matrixes A, B and the vectors X(k), U(k), the following equation (30) is obtained. 
     
       
           X ( k+ 1)= AX ( k )+ BU ( k−d )  (30)  
       
     
     Determining X(k+d) from the equation (30), the following equation (31) is obtained.                X        (     k   +   d     )       =         A   d          X        (   k   )         +     [             A     d   -   1          B             A     d   -   2          B         …         A                 B             B   ]                [           U        (     k   -   1     )                 U        (     k   -   2     )               ⋮             U        (     k   -   d     )             ]                       (   31   )                         
     If matrixes A′ and B′ are defined by the following equations (32), (33), using the model parameters a 1 ′, a 2 ′, b 1 ′, and c 1 ′ which are not subjected to the limit process, a predicted vector XHAT(k+d) is given by the following equation (34).                A   ′     =     [           a1   ′           a2   ′             1       0         ]             (   32   )                 B   ′     =     [           b1   ′           c1   ′             0       0         ]             (   33   )                                       X                 H                 A                   T        (     k   +   d     )         =           A   ′     d          X        (   k   )         +     [             A       ′                 d     -   1            B   ′               A       ′                 d     -   2            B   ′           …           A   ′                     B   ′                 B   ′     ]                [           U        (     k   -   1     )                 U        (     k   -   2     )               ⋮             U        (     k   -   d     )             ]                       (   34   )                         
     The first-row element DTHHAT(k+d) of the predicted vector XHAT(k+d) corresponds to the predicted deviation amount PREDTH(k), and is given by the following equation (35).                      PREDTH        (   k   )       =                DTHHAT        (     k   +   d     )                   =                  α1   ×     DTH        (   k   )         +     α2   ×     DTH        (     k   -   1     )         +                                β                 1   ×     DUT        (     k   -   1     )         +     β                 2   ×     DUT        (     k   -   2     )         +   …   +                                β                 d   ×     DUT        (     k   -   d     )         +     γ                 1     +     γ                 2     +   …   +     γ                 d                     (   35   )                         
     where α 1  represents a first-row, first-column element of the matrix A′ d , α 2  represents a first-row, second-column element of the matrix A′ d , βi represents a first-row, first-column element of the matrix A′ d−i B′, and γi represents a first-row, second-column element of the matrix A′ d−i B′. 
     By applying the predicted deviation amount PREDTH(k) calculated from the equation (35) to the equation (9), and replacing the target values DTHR(k+d+1), DTHR(k+d), and DTHR(k+d−1) respectively with DTHR(k), DTHR(k−1), and DTHR(k−2), the following equation (9a) is obtained. From the equation (9a), the equivalent control input Ueq(k) is calculated.                        D                 U                   T        (   k   )         =                  1   b1          {         (     1   -   a1   -     V                 P                 O                 L                 E       )        P                 R                 E                 D                 T                   H        (   k   )         +                                      (       V                 P                 O                 L                 E     -   a2     )        P                 R                 E                 D                 T                   H        (     k   -   1     )         -   c1   +                                D                 T                 H                   R        (   k   )         +       (       V                 P                 O                 L                 E     -   1     )        D                 T                 H                   R        (     k   -   1     )         -                                V                 P                 O                 L                 E   ×   D                 T                 H                   R        (     k   -   2     )         }                 =                U                 e                   q        (   k   )                            (     9a     )                         
     Using the predicted deviation amount PREDTH(k) calculated from the equation (35), a predicted switching function value σpre(k) is defined by the following equation (36). The reaching law input Urch(k) and the adaptive law input Uadp(k) are calculated respectively from the following equations (10a) and (11a). 
     
       
         σ pre ( k )=( PREDTH ( k )− DTHR ( k− 1))+ VPOLE ( PREDTH ( k− 1)− DTHR ( k− 2))  (36)  
       
       
         
           
             
               
                 
                   
                     U 
                      
                     
                         
                     
                      
                     r 
                      
                     
                         
                     
                      
                     c 
                      
                     
                         
                     
                      
                     
                       h 
                        
                       
                         ( 
                         k 
                         ) 
                       
                     
                   
                   = 
                   
                     
                       
                         - 
                         F 
                       
                       b1 
                     
                      
                     σ 
                      
                     
                         
                     
                      
                     p 
                      
                     
                         
                     
                      
                     r 
                      
                     
                         
                     
                      
                     
                       e 
                        
                       
                         ( 
                         k 
                         ) 
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     10a 
                   
                   ) 
                 
               
             
             
               
                 
                   
                     U 
                      
                     
                         
                     
                      
                     a 
                      
                     
                         
                     
                      
                     d 
                      
                     
                         
                     
                      
                     
                       p 
                        
                       
                         ( 
                         k 
                         ) 
                       
                     
                   
                   = 
                   
                     
                       
                         - 
                         G 
                       
                       b1 
                     
                      
                     
                       
                         ∑ 
                         
                           i 
                           = 
                           0 
                         
                         k 
                       
                        
                       
                         Δ 
                          
                         
                             
                         
                          
                         T 
                          
                         
                             
                         
                          
                         σ 
                          
                         
                             
                         
                          
                         p 
                          
                         
                             
                         
                          
                         r 
                          
                         
                             
                         
                          
                         
                           e 
                            
                           
                             ( 
                             i 
                             ) 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     11a 
                   
                   ) 
                 
               
             
           
         
         
         
             
         
       
     
     The model parameter c 1 ′ is a parameter representing a deviation of the default opening THDEF and disturbance. Therefore, as shown in FIG. 8, the model parameter c 1 ′ changes with disturbance, but can be regarded as substantially constant in a relatively short period. In the present embodiment, the model parameter c 1 ′ is statistically processed, and the central value of its variations is calculated as a default opening deviation thdefadp. The default opening deviation thdefadp is used for calculating the throttle valve opening deviation amount DTH and the target value DTHR. 
     Generally, the method of least squares is known as a method of the statistic process. In the statistic process according to the method of least squares, all data, i.e., all identified parameters c 1 ′, obtained in a certain period are stored in a memory and the stored data is subjected to a batch calculation of the statistic process at a certain timing. However, the batch calculation requires a memory having a large storage capacity for storing all data, and an increased amount of calculations are necessary because inverse matrix calculations are required. 
     Therefore, according to the present embodiment, the sequential method-of-least-squares algorithm for adaptive control which is indicated by the equations (15) through (21) is applied to the statistic process, and the central value of the least squares of the model parameter c 1  is calculated as the default opening deviation thdefadp. 
     Specifically, in the equations (15) through (21), by replacing θ(k) and θ(k) T  with thdefadp, replacing ζ(k) and ζ(k) T  with “1”, replacing ide(k) with ecl(k), replacing KP(k) with KPTH(k), replacing P(k) with PTH(k), and replacing λ 1  and λ 2  respectively with λ 1 ′ and λ 2 ′, the following equations (37) through (40) are obtained. 
           thdefadp ( k+ 1)= thdefadp ( k )+ KPTH ( k ) ec 1( k )  (37)                    K                 P                 T                   H        (   k   )         =       P                 T                   H        (   k   )           1   +     P                 T                   H        (   k   )                     (   38   )                 P                 T                   H        (     k   +   1     )         =       1     λ   1   ′            (     1   -         λ   2   ′        P                 T                   H        (   k   )             λ   1   ′     +       λ   2   ′        P                 T                   H        (   k   )               )        P                 T                   H        (   k   )                 (   39   )                         ec   1 ( k )= c 1′( k )− thdefadp ( k )  ( 40 ) 
     One of the four algorithms described above can be selected according to the setting of the coefficients λ 1 ′ and λ 2 ′. In the equation (39), the coefficient λ 1 ′ is set to a given value other than 0 or 1, and the coefficient λ 2 ′ is set to 1, thus employing the weighted method of least squares. 
     For the calculations of the equations (37) through (40), the values to be stored are thdefadp(k+1) and PTh(k+1) only, and no inverse matrix calculations are required. Therefore, by employing the sequential method-of-least-squares algorithm, the model parameter c 1  can be statistically processed according to the method of least squares while overcoming the shortcomings of a general method of least squares. 
     The default opening deviation thdefadp obtained as a result of the statistic process is applied to the equations (2) and (3), and the throttle valve opening deviation amount DTH(k) and the target value DTHR(k) are calculated from the following equations (41) and (42) instead of the equations (2) and (3). 
     
       
           DTH ( k )= TH ( k )− THDEF+thdefadp   (41)  
       
     
     
       
           DTHR ( k )= THR ( k )− THDEF+thdefadp   (42)  
       
     
     Using the equations (41) and (42), even when the default opening THDEF is shifted from its designed value due to characteristic variations or aging of the hardware, the shift can be compensated to perform an accurate control process. 
     Operation processes executed by the CPU in the ECU  7  for realizing the functions of the adaptive sliding mode controller  21 , the model parameter identifier  22 , and the state predictor  23  will be described below. 
     FIG. 9 is a flowchart showing a process of the throttle valve opening control. The process is executed by the CPU in the ECU  7  in every predetermined period of time (e.g., 2 msec). 
     In step S 11 , a process of setting a state variable shown in FIG. 10 is performed. Calculations of the equations (41) and (42) are executed to determine the throttle valve opening deviation amount DTH(k) and the target value DTHR(k) (steps S 21  and S 22  in FIG.  10 ). The symbol (k) representing a current value may sometimes be omitted as shown in FIG.  10 . 
     In step S 12 , a process of performing calculations of the model parameter identifier as shown in FIG. 11, i.e., a process of calculating the model parameter vector θ(k) from the equation (15a), is carried out. Further, the model parameter vector θ(k) is subjected to the limit process so that the corrected model parameter vector θL(k) is calculated. 
     In step S 13 , a process of performing calculations of the state predictor as shown in FIG. 21 is carried out to calculate the predicted deviation amount PREDTH(k). 
     Next, using the corrected model parameter vector θL(k) calculated in step S 12 , a process of calculating the control input Usl(k) as shown in FIG. 22 is carried out in step S 14 . Specifically, the equivalent control input Ueq, the reaching law input Urch(k), and the adaptive law input Uadp(k) are calculated, and the control input Usl(k) (=duty ratio DUT(k)) is calculated as a sum of these inputs Ueq(k), Urch(k), and Uadp(k). 
     In step S 16 , a process of stability determination of the sliding mode controller as shown in FIG. 29 is carried out. Specifically, the stability is determined based on a differential value of the Lyapunov function, and a stability determination flag FSMCSTAB is set. When the stability determination flag FSMCSTAB is set to “1”, this indicates that the adaptive sliding mode controller  21  is unstable. 
     If the stability determination flag FSMCSTAB is set to “1”, indicating that the adaptive sliding mode controller  21  is unstable, the switching function setting parameter VPOLE is set to a predetermined stabilizing value XPOLESTB (see steps S 231  and S 232  in FIG.  24 ), and the equivalent control input Ueq is set to “0”. That is, the control process by the adaptive sliding mode controller  21  is switched to a control process based on only the reaching law input Urch and the adaptive law input Uadp, to thereby stabilize the control (see steps S 206  and S 208  in FIG.  22 ). 
     Further, when the adaptive sliding mode controller  21  has become unstable, the equations for calculating the reaching law input Urch and the adaptive law input Uadp are changed. Specifically, the values of the reaching law control gain F and the adaptive law control gain G are changed to values for stabilizing the adaptive sliding mode controller  21 , and the reaching law input Urch and the adaptive law input Uadp are calculated without using the model parameter b 1  (see FIGS.  27  and  28 ). According to the above stabilizing process, it is possible to quickly terminate the unstable state of the adaptive sliding mode controller  21 , and to bring the adaptive sliding mode controller  21  back to its stable state. 
     In step S 17 , a process of calculating the default opening deviation thdefadp as shown in FIG. 30 is performed to calculate the default opening deviation thdefadp. 
     FIG. 11 is a flowchart showing a process of performing calculations of the model parameter identifier  22 . 
     In step S 31 , the gain coefficient vector KP(k) is calculated from the equation (20). Then, the estimated throttle valve opening deviation amount DTHHAT(k) is calculated from the equation (18) in step S 32 . In step S 33 , a process of calculating the identifying error idenl(k) as shown in FIG. 12 is carried out. The estimated throttle valve opening deviation amount DTHHAT(k) calculated in step S 32  is applied to the equation (17) to calculate the identifying error ide(k). Further in step S 32 , the dead zone process using the function shown in FIG. 7A is carried out to calculate the corrected identifying error idenl. 
     In step S 34 , the model parameter vector θ(k) is calculated from the equation (15a). Then, the model parameter vector θ(k) is subjected to the stabilization process in step S 35 . That is, each of the model parameters is subjected to the limit process to calculate the corrected model parameter vector θL(k). 
     FIG. 12 is a flowchart showing a process of calculating the identifying error idenl(k) which is carried out in step S 33  shown in FIG.  11 . 
     In step S 51 , the identifying error ide(k) is calculated from the equation (17). Then, it is determined whether or not the value of a counter CNTIDST which is incremented in step S 53  is greater than a predetermined value XCNTIDST that is set according to the dead time d of the controlled object (step S 52 ). The predetermined value XCNTIDST is set, for example, to “3” according to a dead time d=2. Since the counter CNTIDST has an initial value of “0”, the process first goes to step S 53 , in which the counter CNTIDST is incremented by “1”. Then, the identifying error ide(k) is set to “0” in step S 54 , after which the process goes to step S 55 . Immediately after starting identifying the model parameter vector θ(k), no correct identifying error can be obtained by the equation (17). Therefore, the identifying error ide(k) is set to “0” according to steps S 52  through S 54 , instead of using the calculated result of the equation (17). 
     If the answer to the step S 52  is affirmative (YES), then the process immediately proceeds to step S 55 . 
     In step S 55 , the identifying error ide(k) is subjected to a low-pass filtering. Specifically, when identifying the model parameters of an controlled object which has low-pass characteristics, the identifying weight of the method-of-least-squares algorithm for the identifying error ide(k) has frequency characteristics as indicated by the solid line L 1  in FIG.  13 A. By the low-pass filtering of the identifying error ide(k), the frequency characteristics as indicated by the solid line L 1  is changed to a frequency characteristics as indicated by the broken line L 2  where the high-frequency components are attenuated. The reason for executing the low-pass filtering will be described below. 
     The frequency characteristics of the actual controlled object having low-pass characteristics and the controlled object model thereof are represented respectively by the solid lines L 3  and L 4  in FIG.  13 B. Specifically, if the model parameters are identified by the model parameter identifier  22  with respect to the controlled object which has low-pass characteristics (characteristics of attenuating high-frequency components), the identified model parameters are largely affected by the high-frequency-rejection characteristics so that the gain of the controlled object model becomes lower than the actual characteristics in a low-frequency range. As a result, the sliding mode controller  21  excessively corrects the control input. 
     By changing the frequency characteristics of the weighting of the identifying algorithm to the characteristics indicated by the broken line L 2  in FIG. 13A according to the low-pass filtering, the frequency characteristics of the controlled object are changed to frequency characteristics indicated by the broken line L 5  in FIG.  13 B. As a result, the frequency characteristics of the controlled object model is made to coincide with the actual frequency characteristics, or the low frequency gain of the controlled object model is corrected to a level which is slightly higher than the actual gain. Accordinly, it is possible to prevent the control input from being excessively corrected by the sliding mode controller  21 , to thereby improve the robustness of the control system and further stabilize the control system. 
     The low-pass filtering is carried out by storing past values ide(k−i) of the identifying error (e.g., 10 past values for i=1 through 10) in a ring buffer, multiplying the past values by weighting coefficients, and adding the products of the past values and the weighting coefficients. 
     Since the identifying error ide(k) is calculated from the equations (17), (18), and (19), the same effect as described above can be obtained by performing the same low-pass filtering on the throttle valve opening deviation amount DTH(k) and the estimated throttle valve opening deviation amount DTHHAT(k), or by performing the same low-pass filtering on the throttle valve opening deviation amounts DTH(k−1), DTH(k−2) and the duty ratio DUT(k−d−1). 
     Referring back to FIG. 12, the dead zone process as shown in FIG. 14 is carried out in step S 56 . In step S 61  shown in FIG. 14, “n” in the equation (24) is set, for example, to “5” to calculate the square average value DDTHRSQA of an amount of change of the target throttle valve opening THR. Then, an EIDNRLMT table shown in FIG. 15 is retrieved according to the square average value DDTHRSQA to calculate the dead zone width parameter EIDNRLMT (step S 62 ). 
     In step S 63 , it is determined whether or not the identifying error ide(k) is greater than the dead zone width parameter EIDNRLMT. If ide(k) is greater than EIDNRLMT, the corrected identifying error idenl(k) is calculated from the following equation (43) in step S 67 . 
     
       
           idenl ( k )= ide ( k )− EIDNRLMT   (43)  
       
     
     If the answer to step S 63  is negative (NO), it is determined whether or not the identifying error ide(k) is greater than the negative value of the dead zone width parameter EIDNRLMT with a minus sign (step S 64 ). 
     If ide(k) is less than −EIDNRLMT, the corrected identifying error idenl(k) is calculated from the following equation (44) in step S 65 . 
     
       
           idenl ( k )= ide ( k )+ EIDNRLMT   (44)  
       
     
     If the identifying error ide(k) is in the range between +EIDNRLMT and −EIDNRLMT, the corrected identifying error idenl(k) is set to “0” in step S 66 . 
     FIG. 16 is a flowchart showing a process of stabilizing the model parameter vector θ(k), which is carried out in step S 35  shown in FIG.  11 . 
     In step S 71  shown in FIG. 16, flags FA 1 STAB, FA 2 STAB, FB 1 LMT, and FC 1 LMT used in this process are initialized to be set to “0”. In step S 72 , the limit process of the model parameters a 1 ′ and a 2 ′ shown in FIG. 17 is executed. In step S 73 , the limit process of the model parameter b 1 ′ shown in FIG. 19 is executed. In step S 74 , the limit process of the model parameter c 1 ′ shown in FIG. 20 is executed. 
     FIG. 17 is a flowchart showing the limit process of the model parameters a 1 ′ and a 2 ′, which is carried out in the step S 72  shown in FIG.  16 . FIG. 18 is a diagram illustrative of the process shown in FIG. 17, and will be referred to with FIG.  17 . 
     In FIG. 18, combinations of the model parameters a 1 ′ and a 2 ′ which are required to be limited are indicated by “x” symbols, and the range of combinations of the model parameters a 1 ′ and a 2 ′ which are stable are indicated by a hatched region (hereinafter referred to as “stable region”). The limit process shown in FIG. 17 is a process of moving the combinations of the model parameters a 1 ′ and a 2 ′ which are in the outside of the stable region into the stable region at positions indicated by “◯” symbols. 
     In step S 81 , it is determined whether or not the model parameter a 2 ′ is greater than or equal to a predetermined a 2  lower limit value XIDA 2 L. The predetermined a 2  lower limit value XIDA 2 L is set to a negative value greater than “−1”. Stable corrected model parameters a 1  and a 2  are obtained when setting the predetermined a 2  lower limit value XIDA 2 L to “−1”. However, the predetermined a 2  lower limit value XIDA 2 L is set to a negative value greater than “−1” because the matrix A defined by the equation (26) to the “n”th power may occasionally become unstable (which means that the model parameters a 1 ′ and a 2 ′ do not diverge, but oscillate). 
     If a 2 ′ is less than XIDA 2 L in step S 81 , the corrected model parameter a 2  is set to the lower limit value XIDA 2 L, and an a 2  stabilizing flag FA 2 STAB is set to “1”. When the a 2  stabilizing flag FA 2 STAB is set to “1”, this indicates that the corrected model parameter a 2  is set to the lower limit value XIDA 2 L. In FIG. 18, the correction of the model parameter in a limit process P 1  of steps S 81  and S 82  is indicated by the arrow lines with “P 1 ”. 
     If the answer to the step S 81  is affirmative (YES), i.e., if a 2 ′ is greater than or equal to XIDA 2 L, the corrected model parameter a 2  is set to the model parameter a 2 ′ in step S 83 . 
     In steps S 84  and S 85 , it is determined whether or not the model parameter a 1 ′ is in a range defined by a predetermined a 1  lower limit value XIDA 1 L and a predetermined al upper limit value XIDA 1 H. The predetermined al lower limit value XIDA 1 L is set to a value which is equal to or greater than “−2” and lower than “0”, and the predetermined a 1  upper limit value XIDA 1 H is set to “2”, for example. 
     If the answers to steps S 84  and S 85  are affirmative (YES), i.e., if a 1 ′ is greater than or equal to XIDA 1 L and less than or equal to XIDA 1 H, the corrected model parameter a 1  is set to the model parameter a 1 ′ in step S 88 . 
     If a 1 ′ is less than XIDA 1 L in step S 84 , the corrected model parameter a 1  is set to the lower limit value XIDA 1 L and an a 1  stabilizing flag FA 1 STAB is set to “1” in step S 86 . If a 1 ′ is greater than XIDA 1 H in step S 85 , the corrected model parameter a 1  is set to the upper limit value XIDA 1 H and the a 1  stabilizing flag FA 1 STAB is set to “1” in step S 87 . When the a 1  stabilizing flag FA 1 STAB is set to “1”, this indicates that the corrected model parameter a 1  is set to the lower limit value XIDA 1 L or the upper limit value XIDA 1 H. In FIG. 18, the correction of the model parameter in a limit process P 2  of steps S 84  through S 87  is indicated by the arrow lines with “P 2 ”. 
     In step S 90 , it is determined whether or not the sum of the absolute value of the corrected model parameter a 1  and the corrected model parameter a 2  is less than or equal to a predetermined stability determination value XA 2 STAB. The predetermined stability determination value XA 2 STAB is set to a value close to “1” but less than “1” (e.g., “0.99”). 
     Straight lines L 1  and L 2  shown in FIG. 18 satisfy the following equation (45). 
     
       
           a   2 +| a   1 |= XA   2   STAB   (45)  
       
     
     Therefore, in step S 90 , it is determined whether or not the combination of the corrected model parameters a 1  and a 2  is placed at a position on or lower than the straight lines L 1  and L 2  shown in FIG.  18 . If the answer to step S 90  is affirmative (YES), the limit process immediately ends, since the combination of the corrected model parameters a 1  and a 2  is in the stable region shown in FIG.  18 . 
     If the answer to step S 90  is negative (NO), it is determined whether or not the corrected model parameter a 1  is less than or equal to a value obtained by subtracting the predetermined a 2  lower limit value XIDA 2 L from the predetermined stability determination value XA 2 STAB in step S 91  (since XIDA 2 L is less than “0”, XA 2 STAB−XIDA 2 L is greater than XA 2 STAB). If the corrected model parameter a 1  is equal to or less than (XA 2 STAB−XIDA 2 L), the corrected model parameter a 2  is set to (XA 2 STAB−|a 1 |) and the a 2  stabilizing flag FA 2 STAB is set to “1” in step S 92 . 
     If the corrected model parameter a 1  is greater than (XA 2 STAB−XIDA 2 L) in step S 91 , the corrected model parameter a 1  is set to (XA 2 STAB−XIDA 2 L) in step S 93 . Further in step S 93 , the corrected model parameter a 2  is set to the predetermined a 2  lower limit value XIDA 2 L, and the a 1  stabilizing flag FA 1 STAB and the a 2  stabilizing flag FA 2 STAB are set to “1” in step S 93 . 
     In FIG. 18, the correction of the model parameter in a limit process P 3  of steps S 91  and S 92  is indicated by the arrow lines with “P 3 ”, and the correction of the model parameter in a limit processP 4  in steps S 91  and S 93  is indicated by the arrow lines with “P 4 ”. 
     As described above, the limit process shown in FIG. 17 is carried out to bring the model parameters a 1 ′ and a 2 ′ into the stable region shown in FIG. 18, thus calculating the corrected model parameters a 1  and a 2 . 
     FIG. 19 is a flowchart showing a limit process of the model parameter b 1 ′, which is carried out in step S 73  shown in FIG.  16 . 
     In steps S 101  and S 102 , it is determined whether or not the model parameter b 1 ′ is in a range defined by a predetermined b 1  lower limit value XIDB 1 L and a predetermined b 1  upper limit value XIDB 1 H. The predetermined b 1  lower limit value XIDB 1 L is set to a positive value (e.g., “0.1”), and the predetermined b 1  upper limit value XIDB 1 H is set to “1”, for example. 
     If the answers to steps S 101  and S 102  are affirmative (YES), i.e., if b 1 ′ is greater than or equal to XIDB 1 L and less than or equal to XIDB 1 H, the corrected model parameter b 1  is set to the model parameter b 1 ′ in step S 105 . 
     If b 1 ′ is less than XIDB 1 L in step S 101 , the corrected model parameter b 1  is set to the lower limit value XIDB 1 L, and a b 1  limiting flag FB 1 LMT is set to “1” in step S 104 . If b 1 ′ is greater than XIDB 1 H in step S 102 , then the corrected model parameter b 1  is set to the upper limit value XIDB 1 H, and the b 1  limiting flag FB 1 LMT is set to “1” in step S 103 . When the b 1  limiting flag FB 1 LMT is set to “1”, this indicates that the corrected model parameter b 1  is set to the lower limit value XIDB 1 L or the upper limit value XIDB 1 H. 
     FIG. 20 is a flowchart showing a limit process of the model parameter c 1 ′, which is carried out in step S 74  shown in FIG.  16 . 
     In steps S 111  and S 112 , it is determined whether or not the model parameters c 1 ′ is in a range defined by a predetermined c 1  lower limit value XIDC 1 L and a predetermined c 1  upper limit value XIDC 1 H. The predetermined c 1  lower limit value XIDC 1 L is set to “−60”, for example, and the predetermined c 1  upper limit value XIDC 1 H is set to “60”, for example. 
     If the answers to steps S 111  and S 112  are affirmative (YES), i.e., if c 1 ′ is greater than or equal to XIDC 1 L and less than or equal to XIDC 1 H, the corrected model parameter c 1  is set to the model parameter c 1 ′ in step S 115 . 
     If c 1 ′ is less than XIDC 1 L in step S 111 , the corrected model parameter c 1  is set to the lower limit value XIDC 1 L, and a c 1  limiting flag FC 1 LMT is set to “1” in step S 114 . If c 1 ′ is greater than XIDC 1 H in step S 112 , the corrected model parameter c 1  is set to the upper limit value XIDC 1 H, and the c 1  limiting flag FC 1 LMT is set to “1” in step S 113 . When the c 1  limiting flag FC 1 LMT is set to “1”, this indicates that the corrected model parameter c 1  is set to the lower limit value XIDC 1 L or the upper limit value XIDC 1 H. 
     FIG. 21 is a flowchart showing a process of calculations of the state predictor, which is carried out in step S 13  shown in FIG.  9 . 
     In step S 121 , the matrix calculations are executed to calculate the matrix elements α 1 , α 2 , β 1 , β 2 , γ 1  through γd in the equation (35). 
     In step S 122 , the predicted deviation amount PREDTH(k) is calculated from the equation (35). 
     FIG. 22 is a flowchart showing a process of calculation of the control input Usl (=DUT) applied to the throttle valve actuating device  10 , which is carried out in step S 14  shown in FIG.  9 . 
     In step S 201 , a process of calculation of the predicted switching function value σpre, which is shown in FIG. 23, is executed. In step S 202 , a process of calculation of the integrated value of the predicted switching function value σpre, which is shown in FIG. 26, is executed. In step S 203 , the equivalent control input Ueq is calculated from the equation (9). In step S 204 , a process of calculation of the reaching law input Urch, which is shown in FIG. 27, is executed. In step S 205 , a process of calculation of the adaptive law input Uadp, which is shown in FIG. 28, is executed. 
     In step S 206 , it is determined whether or not the stability determination flag FSMCSTAB set in the process shown in FIG. 29 is “1”. When the stability determination flag FSMCSTAB is set to “1”, this indicates that the adaptive sliding mode controller  21  is unstable. 
     If FSMCSTAB is “0” in step S 206 , indicating that the adaptive sliding mode controller  21  is stable, the control inputs Ueq, Urch, and Uadp which are calculated in steps S 203  through S 205  are added, thereby calculating the control input Usl in step S 207 . 
     If FSMCSTAB is “1” in step S 206 , indicating that the adaptive sliding mode controller  21  is unstable, the sum of the reaching law input Urch and the adaptive law input Uadp is calculated as the control input Usl. In other words, the equivalent control input Ueq is not used for calculating the control input Usl, thus preventing the control system from becoming unstable. 
     In steps S 209  and S 210 , it is determined whether or not the calculated control input Usl is in a range defined between a predetermined upper limit value XUSLH and a predetermined lower limit value XUSLL. If the control input Usl is in the range between XUSLH and XUSLL, the process immediately ends. If the control input Usl is less than or equal to the predetermined lower limit value XUSLL in step S 209 , the control input Usl is set to the predetermined lower limit value XUSLL in step S 212 . If the control input Usl is greater than or equal to the predetermined upper limit value XUSLH in step S 210 , the control input Usl is set to the predetermined upper limit value XUSLH in step S 211 . 
     FIG. 23 is a flowchart showing a process of calculating the predicted switching function value σpre, which is carried out in step S 201  shown in FIG.  22 . 
     In step S 221 , the process of calculating the switching function setting parameter VPOLE, which is shown in FIG. 24 is executed. Then, the predicted switching function value σpre(k) is calculated from the equation (36) in step S 222 . 
     In steps S 223  and S 224 , it is determined whether or not the calculated predicted switching function value σpre(k) is in a range defined between a predetermined upper value XSGMH and a predetermined lower limit value XSGML. If the calculated predicted switching function value σpre(k) is in the range between XSGMH and XSGML, the process shown in FIG. 23 immediately ends. If the calculated predicted switching function value σpre(k) is less than or equal to the predetermined lower limit value XSGML in step S 223 , the calculated predicted switching function value σpre(k) is set to the predetermined lower limit value XSGML in step S 225 . If the calculated predicted switching function value σpre(k) is greater than or equal to the predetermined upper limit value XSGMH in step S 224 , the calculated predicted switching function value σpre(k) is set to the predetermined upper limit value XSGMH in step S 226 . 
     FIG. 24 is a flowchart showing a process of calculating the switching function setting parameter VPOLE, which is carried out in step S 221  shown in FIG.  23 . 
     In step S 231 , it is determined whether or not the stability determination flag FSMCSTAB is “1”. If FSMCSTAB is “1” in step S 231 , indicating that the adaptive sliding mode controller  21  is unstable, the switching function setting parameter VPOLE is set to a predetermined stabilizing value XPOLESTB in step S 232 . The predetermined stabilizing value XPOLESTB is set to a value which is greater than “−1” but very close to “−1” (e.g., “−0.999”). 
     If FSMCSTAB is “0”, indicating that the adaptive sliding mode controller  21  is stable, an amount of change DDTHR(k) in the target value DTHR(k) is calculated from the following equation (46) in step S 233 . 
     
       
           DDTHR ( k )= DTHR ( k )− DTHR ( k− 1)  (46)  
       
     
     In step S 234 , a VPOLE map is retrieved according to the throttle valve opening deviation amount DTH and the amount of change DDTHR calculated in step S 233  to calculate the switching function setting parameter VPOLE. As shown in FIG. 25A, the VPOLE map is set so that the switching function setting parameter VPOLE increases when the throttle valve opening deviation amount DTH has a value in the vicinity of “0”, i.e., when the throttle valve opening TH is in the vicinity of the default opening THDEF, and the switching function setting parameter VPOLE has a substantially constant value regardless of changes of the throttle valve opening deviation amount DTH, when the throttle valve opening deviation amount DTH has values which are not in the vicinity of “0”. The VPOLE map is also set so that the switching function setting parameter VPOLE increases as the amount of change DDTHR in target value increases as indicated by the solid line in FIG. 25B, and the switching function setting parameter VPOLE increases as the amount of change DDTHR in the target value has a value in the vicinity of “0” as indicated by the broken line in FIG. 25B, when the throttle valve opening deviation amount DTH has a value in the vicinity of “0”. 
     Specifically, when the target value DTHR for the throttle valve opening changes greatly in the decreasing direction, the switching function setting parameter VPOLE is set to a relatively small value. This makes it possible to prevent the throttle valve  3  from colliding with the stopper for stopping the throttle valve  3  in the fully closed position. In the vicinity of the default opening THDEF, the switching function setting parameter VPOLE is set to a relatively large value, which improves the controllability in the vicinity of the default opening THDEF. 
     As shown in FIG. 25C, the VPOLE map may be set so that the switching function setting parameter VPOLE decreases when the throttle valve opening TH is in the vicinity of the fully closed opening or the fully open opening. Therefore, when the throttle valve opening TH is in the vicinity of the fully closed opening or the fully open opening, the speed for the throttle valve opening TH to follow up the target opening THR is reduced. As a result, collision of the throttle valve  3  with the stopper can more positively be avoided (the stopper also stops the throttle valve  3  in the fully open position). 
     In steps S 235  and S 236 , it is determined whether or not the calculated switching function setting parameter VPOLE is in a range defined between a predetermined upper limit value XPOLEH and a predetermined lower limit XPOLEL. If the switching function setting parameter VPOLE is in the range between XPOLEH and XPOLEL, the process shown immediately ends. If the switching function setting parameter VPOLE is less than or equal to the predetermined lower limit value XPOLEL in step S 236 , the switching function setting parameter VPOLE is set to the predetermined lower limit value XPOLEL in step S 238 . If the switching function setting parameter VPOLE is greater than or equal to the predetermined upper limit value XPOLEH in step S 235 , the switching function setting parameter VPOLE is set to the predetermined upper limit value XPOLEH in step S 237 . 
     FIG. 26 is a flowchart showing a process of calculating an integrated value of σpre, SUMSIGMA, of the predicted switching function value σpre. This process is carried out in step S 202  shown in FIG.  22 . The integrated value SUMSIGMA is used for calculating the adaptive law input Uadp in the process shown in FIG. 28 which will be described later (see the equation (11a)). 
     In step S 241 , the integrated value SUMSIGMA is calculated from the following equation (47) where ΔT represents a calculation period. 
     
       
           SUMSIGMA ( k )= SUMSIGMA ( k− 1)+σ pre×ΔT   (47)  
       
     
     In steps S 242  and S 243 , it is determined whether or not the calculated integrated value SUMSIGMA is in a range defined between a predetermined upper limit value XSUMSH and a predetermined lower limit value XSUMSL. If the integrated value SUMSIGMA is in the range between XSUMSH and XSUMSL, the process immediately ends. If the integrated value SUMSIGMA is less than or equal to the predetermined lower limit value XSUMSL in step S 242 , the integrated value SUMSIGMA is set to the predetermined lower limit value XSUMSL in step S 244 . If the integrated value SUMSIGMA is greater than or equal to the predetermined upper limit value XSUMSH in step S 243 , the integrated value SUMSIGMA is set to the predetermined upper limit value XSUMSH in step S 245 . 
     FIG. 27 is a flowchart showing a process of calculating the reaching law input Urch, which is carried out in step S 204  shown in FIG.  22 . 
     In step S 261 , it is determined whether or not the stability determination flag FSMCSTAB is “1”. If the stability determination flag FSMCSTAB is “0”, indicating that the adaptive sliding mode controller  21  is stable, the control gain F is set to a normal gain XKRCH in step S 262 , and the reaching law input Urch is calculated from the following equation (48), which is the same as the equation (10a), in step S 263 . 
     
       
           Urch=−F×σpre/b   1   (48)  
       
     
     If the stability determination flag FSMCSTAB is “1”, indicating that the adaptive sliding mode controller  21  is unstable, the control gain F is set to a predetermined stabilizing gain XKRCHSTB in step S 264 , and the reaching law input Urch is calculated according to the following equation (49), which does not include the model parameter b 1 , in step S 265 . 
     
       
           Urch=−F×σpre   (49)  
       
     
     In steps S 266  and S 267 , it is determined whether the calculated reaching law input Urch is in a range defined between a predetermined upper limit value XURCHH and a predetermined lower limit value XURCHL. If the reaching law input Urch is in the range between XURCHH and XURCHL, the process immediately ends. If the reaching law input Urch is less than or equal to the predetermined lower limit value XURCHL in step S 266 , the reaching law input Urch is set to the predetermined lower limit value XURCHL in step S 268 . If the reaching law input Urch is greater than or equal to the predetermined upper limit value XURCHH in step S 267 , the reaching law input Urch is set to the predetermined upper limit value XURCHH in step S 269 . 
     As described above, when the adaptive sliding mode controller  21  becomes unstable, the control gain F is set to the predetermined stabilizing gain XKRCHSTB, and the reaching law input Urch is calculated without using the model parameter b 1 , which brings the adaptive sliding mode controller  21  back to its stable state. When the identifying process carried out by the model parameter identifier  22  becomes unstable, the adaptive sliding mode controller  21  becomes unstable. Therefore, by using the equation (49) that does not include the model parameter b 1  which has become unstable, the adaptive sliding mode controller  21  can be stabilized. 
     FIG. 28 is a flowchart showing a process of calculating the adaptive law input Uadp, which is carried out in step S 205  shown in FIG.  22 . 
     In step S 271 , it is determined whether or not the stability determination flag FSMCSTAB is “1”. If the stability determination flag FSMCSTAB is “0”, indicating that the adaptive sliding mode controller  21  is stable, the control gain G is set to a normal gain XKADP in step S 272 , and the adaptive law input Uadp is calculated from the following equation (50), which corresponds to the equation (11a), in step S 273 . 
     
       
           Uadp=−G×SUMSIGMA/b   1   (50)  
       
     
     If the stability determination flag FSMCSTAB is “1”, indicating that the adaptive sliding mode controller  21  is unstable, the control gain G is set to a predetermined stabilizing gain XKADPSTB in step S 274 , and the adaptive law input Uadp is calculated according to the following equation (51), which does not include the model parameter b 1 , in step S 275 . 
     
       
           Uadp=−G×SUMSIGMA   (51)  
       
     
     As described above, when the adaptive sliding mode controller  21  becomes unstable, the control gain G is set to the predetermined stabilizing gain XKADPSTB, and the adaptive law input Uadp is calculated without using the model parameter b 1 , which brings the adaptive sliding mode controller  21  back to its stable state. 
     FIG. 29 is a flowchart showing a process of determining the stability of the sliding mode controller, which is carried out in step S 16  shown in FIG.  9 . In this process, the stability is determined based on a differential value of the Lyapunov function, and the stability determination flag FSMCSTAB is set according to the result of the stability determination. 
     In step S 281 , a switching function change amount Dσpre is calculated from the following equation (52). Then, a stability determining parameter SGMSTAB is calculated from the following equation (53) in step S 282 . 
     
       
           Dσpre=σpre ( k )−σ pre ( k− 1)  (52)  
       
     
     
       
           SGMSTAB=Dσpre×σpre ( k )  (53)  
       
     
     In step S 283 , it is determined whether or not the stability determination parameter SGMSTAB is less than or equal to a stability determining threshold XSGMSTAB. If SGMSTAB is greater than XSGMSTAB, it is determined that the adaptive sliding mode controller  21  may possibly be unstable, and an unstability detecting counter CNTSMCST is incremented by “1” in step S 285 . If SGMSTAB is less than or equal to XSGMSTAB, the adaptive sliding mode controller  21  is determined to be stable, and the count of the unstability detecting counter CNTSMCST is not incremented but maintained in step S 284 . 
     In step S 286 , it is determined whether or not the value of the unstability detecting counter CNTSMCST is less than or equal to a predetermined count XSSTAB. If CNTSMCST is less than or equal to XSSTAB, the adaptive sliding mode controller  21  is determined to be stable, and a first determination flag FSMCSTAB 1  is set to “0” in step S 287 . If CNTSMCST is greater than XSSTAB, the adaptive sliding mode controller  21  is determined to be unstable, and the first determination flag FSMCSTAB 1  is set to “1” in step S 288 . The value of the unstability detecting counter CNTSMCST is initialized to “0”, when the ignition switch is turned on. 
     In step S 289 , a stability determining period counter CNTJUDST is decremented by “1”. It is determined whether or not the value of the stability determining period counter CNTJUDST is “0” in step S 290 . The value of the stability determining period counter CNTJUDST is initialized to a predetermined determining count XCJUDST, when the ignition switch is turned on. Initially, therefore, the answer to step S 290  is negative (NO), and the process immediately goes to step S 295 . 
     If the count of the stability determining period counter CNTJUDST subsequently becomes “0”, the process goes from step S 290  to step S 291 , in which it is determined whether or not the first determination flag FSMCSTAB 1  is “1”. If the first determination flag FSMCSTAB 1  is “0”, a second determination flag FSMCSTAB 2  is set to “0” in step S 293 . If the first determination flag FSMCSTAB 1  is “1”, the second determination flag FSMCSTAB 2  is set to “1” in step S 292 . 
     In step S 294 , the value of the stability determining period counter CNTJUDST is set to the predetermined determining count XCJUDST, and the unstability detecting counter CNTSMCST is set to “0”. Thereafter, the process goes to step S 295 . 
     In step S 295 , the stability determination flag FSMCSTAB is set to the logical sum of the first determination flag FSMCSTAB 1  and the second determination flag FSMCSTAB 2 . The second determination flag FSMCSTAB 2  is maintained at “1” until the value of the stability determining period counter CNTJUDST becomes “0”, even if the answer to step S 286  becomes affirmative (YES) and the first determination flag FSMCSTAB 1  is set to “0”. Therefore, the stability determination flag FSMCSTAB is also maintained at “1” until the value of the stability determining period counter CNTJUDST becomes “0”. 
     FIG. 30 is a flowchart showing a process of calculating the default opening deviation thdefadp, which is carried out in step S 17  shown in FIG.  9 . 
     In step S 251  shown FIG. 30, a gain coefficient KPTH(k) is calculated according to the following equation (54). 
     
       
           KPTH ( k )= PTH ( k− 1)/(1 +PTH ( k− 1))  (54)  
       
     
     where PTH(k−1) represents a gain parameter calculated in step S 253  when the present process was carried out in the preceding cycle. 
     In step S 252 , the model parameter c 1 ′ calculated in the process of calculations of the model parameter identifier as shown in FIG.  11  and the gain coefficient KPTH(k) calculated in step S 251  are applied to the following equation (55) to calculate a default opening deviation thdefadp(k). 
     
       
           thdefadp ( k )= thdefadp ( k− 1)+ KPTH ( k )×( c   1 ′− thdefadp ( k− 1))  (55)  
       
     
     In step S 253 , a gain parameter PTH(k) is calculated from the following equation (56): 
     
       
           PTH ( k )={1 −PTH ( k− 1)}/( XDEFADPW+PTH ( k− 1))}× PTH ( k− 1)/ XDEFADPW   (56)  
       
     
     The equation (56) is obtained by setting λ 1 ′ and λ 2 ′ in the equation (39) respectively to a predetermined value XDEFADP and “1”. 
     According to the process shown in FIG. 30, the model parameter c 1 ′ is statistically processed by the sequential method-of-weighted-least-squares algorithm to calculate the default opening deviation thdefadp. 
     In the present embodiment, the model parameter identifier  22  shown in FIG. 3 corresponds to an identifying means, and the state predictor  23  shown in FIG. 3 corresponds to a predicting means. 
     Although a certain preferred embodiment of the present invention has been shown and described in detail, it should be understood that various changes and modifications may be made therein without departing from the scope of the appended claims.