Abstract:
A plant model is constructed with a mathematical procedure known as the Kalman filter. The presumed value or predictive value of a certain parameter of a plant is compared with the actually-measured value of the parameter, and the state transition matrix and drive matrix of the Kalman filter are corrected. When a parameter which is a direct controlled object is difficult of actual measurement, the model is corrected by utilizing another parameter whose actual measurement value is easily obtained. When the plant has a feedback control system based on the proportional-plus-integral control, the gain of the system is corrected in correspondence with the corrections of the elements of the matrices.

Description:
This is a continuation of application Ser. No. 280,206 filed July 6, 1981, now abandoned. 
    
    
     BACKGROUND OF THE INVENTION 
     This invention relates to a plant control method. 
     In order to operate, for example, a thermal power plant stably and at high efficiency, it is important that the temperature of main steam generated by a boiler is kept constant without being affected even by the disturbances of a load change etc. In starting the plant, the main steam temperature needs to be raised correctly in accordance with a predetermined heating curve. 
     To the end of enhancing the controllability, various systems have heretofore been developed. The inventors of the present application have also proposed the system of U.S. Pat. No. 3,877,636. The patented invention derives the desired value and predictive value of the main steam temperature after a predetermined time by means of a computer which mathematically models the the characteristics of a boiler, and controls the starting of the boiler so as to render the difference of both the desired and predictive values null. Although this system has achieved a certain result in the enhancement of the controllability, a more enhanced performance is desired. Further, the proposed system needs to be made applicable also to the steady load operation. 
     One of the reasons why a fully satisfactory performance is not obtained even with the aforecited invention is that the characteristics of the boiler cannot be correctly simulated. Another reason is that the characteristics vary depending upon the operating conditions. 
     SUMMARY OF THE INVENTION 
     An object of this invention is to propose a plant control method in which a model for simulating the characteristics of a plant is corrected with the actually-measured values of the plant in an on-line fashion, thereby to reproduce the plant characteristics more faithfully and to obtain more appropriate predictive values. 
     Another object of this invention is to provide a specific measure for realizing the plant control method of this invention for a boiler plant, the measure employing a secondary superheater as a model representative of the plant. 
     Still another object of this invention is to provide a specific measure for realizing the plant control method of this invention for a boiler plant whose fuel is coal, the measure employing a model of a coal mill conjointly. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a schematic flow diagram showing a conceptual construction in which this invention is applied to a coal-burning thermal power plant. 
     FIG. 2 is a block diagram for explaining the fundamental concept of this invention according to which the model of a plant is constructed and corrected by exploiting the idea of a Kolman filter. 
     FIG. 3 is a graph for explaining the idea of the model correction in FIG. 2. 
     FIG. 4 is a simplified structural view of a secondary superheater. 
     FIG. 5 is a simplified structural view of a mill system. 
     FIG. 6 is a flow diagram showing control blocks in the application of this invention to a thermal power plant. 
     FIG. 7 is a waveform diagram for explaining the effect of this invention. 
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT 
     FIG. 1 shows a conceptual construction at the time when this invention is applied to a coal-burning thermal power plant. Numeral 10 designates a console, numeral 20 a computer, numeral 30 a coal mill system, numeral 40 a boiler system, and numeral 50 a turbine generator system. 
     The operator of the thermal power plant executes necessary operations on the console 10 with reference to data in various parts of the plant transmitted through the computer 20, and data transmitted from a higher-ranked control station such as central load-dispatching office which is not shown, etc. The computer 20 gives the various parts of the plant necessary control signals which have been prepared from the data of the various parts of the plant, signals responsive to the operations of the console 10, etc. in accordance with a program stored in advance. The coal mill system 30 is constructed of a coal bunker 301, a coal feed belt 302, a belt driving motor 303, a mill 310, blowers 321 and 322, and dampers 323 and 324. The mill 310 consists of a table 311, a table driving motor 312, and a plurality of balls 313. Coal supplied from the bunker 301 through the belt 302 into the mill 310 is pulverized between the table 311 and the balls 313. Air is blown into the mill 310 by the blowers 321 and 322, and the coal pulverized into fine powdery coal is conveyed by the air to the burner of a boiler to be described later. In order to control the quantity of the coal, the motor 303 is controlled by the computer 20. In order to control the quantity of the air, the dampers 323 and 324 are controlled by respective actuators 331 and 332 which are operated by the outputs of the computer 20. An air pipe is branched between the blower 321 and the damper 324, to feed air which is supplied around the burner of the boiler. Data relevant to the mill, the quantity of the air to be blown into the mill (F MA ), a mill differential pressure (ΔP M ) being the pressure difference between a mill inlet pressure and a mill internal pressure, and the quantity of the air to be fed around the burner (F BA ) are respectively detected by sensors 341, 342 and 343, and the detected signals are transmitted to the computer 20. 
     The boiler system 40 is constructed of a feed water pump 401, a feed water control valve 402, an evaporator 403, a primary superheater 404, a secondary superheater 405, a gas recirculating blower 406 and a burner 407. Water supplied from the feed water pump 401 is turned into steam by the evaporator 403. The steam is turned by the primary and secondary superheaters 404 and 405 into main steam which is superheated steam and which is supplied to the turbine generator system 50 to be described later. The quantity of heat supplied by the burner 407 is used for turning the water into the steam by the evaporator and turning the steam into the superheated steam by the superheaters, and partially escapes from a chimney into the atmospheric air along with an emission gas. Part of the gas to be emitted from the chimney is returned to the boiler by the gas recirculating blower 406. In order to control the flow rate of the steam which the boiler supplies, the control valve 402 is operated by the output of the computer 20. As data relevant to the boiler for controlling the temperature of the main steam, the quantity of the feed water, the outlet temperature of the primary superheater, the outlet pressure thereof, the outlet steam flow rate thereof, the main steam temperature and the gas recirculation flow rate are respectively detected by sensors 411, 412, 413, 414, 415 and 416, and the detected signals are transmitted to the computer 20. 
     The turbine generator system 50 is constructed of a turbine control valve 501, a high-pressure turbine 502, a medium- and low-pressure turbine 503, a steam condenser 504, and a generator 505 directly coupled to a turbine rotor. The control valve 501 is operated by the output of the computer 20, and the main steam in an amount corresponding to the degree of opening of the valve is supplied to the turbines 502 and 503 to rotate the generator 505. The steam having finished its work in the turbines is turned back into water by the condenser 504. This water is supplied to the boiler again by the feed water pump 401 stated before. A generator output is detected by a sensor 509, and the detected signal is transmitted to the computer 20. 
     The computer 20 is given various requests concerning the plant operation from the console 10, and delivers predetermined control signals on the basis of the data received from the plant and the program stored in advance. In the embodiment of FIG. 1, the computer 20 contains therein a computation model of the mill system 30 and a computation model of the boiler system 40. Each of the models utilizes a technique known as the Kalman filter, and has its precision enhanced. In the present embodiment, the models of the mill system 30 and the boiler system 40 are formed as being independent as will be described later. As to the mill system 30, the model is corrected on the basis of the actually-measured values and presumed values of the feed air quantity and the mill differential pressure, and the quantity of the fine powdery coal to be supplied to the boiler burner 407 is presumed. As to the boiler system 40, the model is corrected on the basis of the actually-measured value and presumed value of the main steam temperature, and the predictive value of the main steam temperature after a predetermined time is derived. 
     The model itself of this invention and the construction of a control system employing the model will be specifically described later. 
     FIG. 2 is a block diagram for explaining the fundamental concept according to which the model of the plant is constructed and corrected by exploiting the idea of the Kalman filter in this invention. 
     Numeral 201 designates the controlled plant or system for which the model is to be constructed, and the dynamic characteristic of which is supposed to be expressed by the following state transition equation: 
     
          (i+1)=φ(i)· (i)+H(i)· (i)           (1) 
    
     where 
      (i): value of an n-dimensional state vector at a time i. 
     φ(i): (n×n)-dimensional state transition matrix. 
     H(i): (n×r)-dimensional drive matrix. 
      (i): r-dimensional probability variable vector representative of system noise. 
     The noise vector  (i) is assumed to be the white random noise which is its average value and dispersion expressed as follows: 
     
         E{ (i)}= (i)                                               (2) 
    
     
         E{ (i)- (i)}{ (j)- (j)}&#39;=δ.sub.ij U                  (3) 
    
     where 
     U: (r×r)-dimensional positive fixed value matrix. ##EQU1## {}&#39;: transposition. Numerals 203 and 205 indicate measuring instruments, which derive the state vector  (i) and the noise vector  (i) in Equation (1) respectively. Here, the measuring instrument 203 is supposed to have the following observation equation: 
     
          (i)=C(i)· (i)+ (i)                               (4) 
    
     where 
      (i): m-dimensional observation vector. 
     C(i): (m×n)-dimensional observation matrix. 
      (i): m-dimensional observation noise vector. 
     The observation noise vector  (i) is also supposed to be the white random noise and to have properties similar to those indicated in Equations (2) and (3). In addition, the observation noise vector  (i) is assumed to be independent of the system noise vector  (i) and the initial value vector (O). Stated differently, the noise vector  (i) is a kind of disturbance acting on the system. 
     Numeral 207 indicates a numerical expression model of the plant, the state transition equation of which is given by the aforecited equation (1). It is the so-called Kalman filter that executes the following mathematical operation by exerting the most probable presumed value   of the state vector  (i) on the state transition equation φ(i) and exerting the noise vector  (i) on the drive matrix H(i): ##EQU2## where   denotes an error vector, and  (i) the presumed quantity of the model. ##EQU3## Numerals 209 and 211 in FIG. 2 indicate an error variance and an observation matrix for calculating Equation (5), respectively. C(i) of the observation matrix 211 is identical to the observation matrix C(i) of the measuring instrument 203. 
     Numeral 213 indicates a numerical expression model of the plant. Although this model is essentially the same as the model 207, it delivers the most probable predictive value  (i, j) after j samplings instead of the most probable presumed value  (i) at the present time (i). That is: 
     
          (i, j)= (i, j)+ (i)                                       (10) 
    
     where 
     
          (i, j)=φ(i) (i, j-1)+H(i) (i)                         (11) 
    
     j=1, 2, . . . n 
     The other conditions are the same as in the aforecited operation as the Kalman filter. 
     A block 215 indicated by a broken line in FIG. 2 is another Kalman filter. It is meant that a most probable presumed value by the Kalman filter 215 may well be used as one of noise vectors in the Kalman filters indicated blocks 207 and 213. In order to simplify the state transition equation of the plant, the plant may be considered a combination of appropriate units into which the plant is divided. In the example of FIG. 1 in which this invention is applied to the plant, the boiler system is made the model unit for predicting the main steam temperature, while the mill system is made the model unit for presuming the coal quantity to be supplied from the mill to the boiler. The details of this fact will be described later. 
     Numeral 217 indicates a model correcting function provided in this invention. It serves to continually correct the model in order to attain a higher control precision when the plant is controlled by utilizing the mathematical method known as the Kalman filter. An example of the model correcting function will be described with reference to FIG. 3. 
     In FIG. 3, the abscissa represents the time, and t i , t i+1  and t i+2  indicate the sampling times of the state vector   as are respectively spaced time intervals Δt. The ordinate represents the actually-measured value and predictive value of the state vector. As given by Equation (1), the state transition equation of the plant is as follows: 
     
          (i+1)=φ(i) (i)+H(i) (i) 
    
     At the present time t i , the predictive values at the times t i+1  and t i+2  or vectors   p  (i, 1) and   p  (i, 2) become as follows respectively: 
     
         .sub.p (i, 1)=φ(i) (i)+H(i) (i)                       (12) 
    
     
         .sub.p (i, 2)=φ(i) (i, 1)+H(i) (i)                    (13) 
    
     On the other hand, letting   a  (i+1) and   a  (i+2) denote the actual observation values from the plant at the respective times, the corresponding error vectors   are as follows: 
     
          (i+1)= .sub.a (i+1)- .sub.p (i, 1)                        (14) 
    
     
          (i+2)= .sub.a (i+2)- .sub.p (i, 2)                        (15) 
    
     By subtracting (14) from (15) and then substituting (12) and (13), ##EQU4## If the state transition equation φ(i) is the optimum, { (i+1)- (i+2)} in Equation (17) ought to have become zero. Letting φ MS  denote a state transition matrix at this time, ##EQU5## The error vector φ&#39;(i+1) at the time when the prediction is made by the use of the matrix φ MS  becomes as follows: 
     
         e.sub.MS (i+1)= .sub.a (i+1)-φ.sub.MS (i) (i)-H(i) (i) (19) 
    
     If the drive matrix H(i) is the optimum, the error vector  (i+1) ought to have become zero. Letting H MS  (i) denote the drive matrix at this time, ##EQU6## 
     That is, every sampling of the state vector, the actually-measured value and the predictive value are compared, and the state transition matrix φ MS  (i) and the drive matrix H MS  (i) respectively satisfying Equation (18) and Equation (20) are introduced into the models 207 and 213, whereby the most probable presumed value  (i) and most probable predictive value  (i, j) at higher precision can be obtained. 
     In order to solve Equations (18) and (20), repeated operations are necessary. Accordingly, not only the computer needs to be large-sized, but also a long period of time is taken for the operations. In order to avoid this drawback, that one of the elements of the system which is the most liable to involve errors in modeling may be noticed and corrected. When, for example, the boiler system is considered, the heat transfer rate α in the state transition matrix φ(i) and the quantity of heat generation H u  in the drive matrix H(i) shall be respectively corrected by integral calculations as follows: 
     
         α=∫{ (i+2)- (i+1)}dt                            (21) 
    
     
         H.sub.u =∫ .sub.MS (i+1)dt                            (22) 
    
     While the above model corrections are based on the predictive values, presumed values may well be relied on as stated below. 
     The error vector is calculated from the presumed value and the actually-measured value every sampling, and it is handled similarly to the error vector in Equations (21) and (22) so as to correct the heat transfer rate α and the quantity of heat generation H u . The model which requires to prediction calculation is preferably constructed in this manner because the operation time for the prediction can be saved. The presumed values and error vectors in this case are as follows, and can be processed in the same way as in the foregoing idea: 
     
         .sub.m (i-1)=φ(i-2) (i-2)+H(i-2) (i-2) 
    
     
         .sub.m (i)=φ(i-2) (i-1)+H(i-2) (i-2) 
    
     
          (i-1)= .sub.a (i-1)- .sub.m (i-1) 
    
     
          (i)= .sub.a (i)- .sub.m (i) 
    
     There will now be described specific examples in which the fundamental idea of this invention above described is applied to the thermal power plant shown in FIG. 1. First, the application to the prediction of the main steam temperature in the boiler system will be described. 
     FIG. 4 is a simplified structural view of the secondary superheater in which the secondary superheater is considered as a concentrated constant in order to simplify the model. In the figure, numeral 421 designates the metal pipe of the secondary superheater, and numeral 422 an outer wall. The main steam flows inside the metal pipe, and the gas flows between the pipe and the outer wall. Here, the following equation holds: 
     From the law of conservation of energy, ##EQU7## where V S2  SH ; volume of the internal fluid of the secondary superheater. 
     γ S2  SH ; specific weight of the internal fluid of the secondary superheater. 
     H S2  SH ; enthalpy of the internal fluid of the secondary superheater. 
     F S2  SH ; flow rate of the internal fluid of the secondary superheater. 
     θ S2  SH ; temperature of the internal fluid at the outlet of the secondary superheater (steam temperature). 
     θ m2  SH ; average temperature of the metal of the secondary superheater. 
     A S2  SH ; heat transfer area of the secondary superheater. 
     α mS2SH  ; rate of heat transfer from the secondary superheater to the internal fluid 
     
         =α.sub.mS2SH, r.sup.(F S2SH/F.sub.S2SH, r).sup.0.8   (24) 
    
     α mS2SH , r : rate of heat transfer from the metal to the internal fluid in a rated state. 
     F  S2SH , r ; flow rate of the internal fluid of the secondary superheater in the rated state. 
     H  S1SH  ; enthalpy of the internal fluid at the inlet of the secondary superheater. 
     On the other hand, for the metal of the secondary superheater, the following equation holds from the law of conservation of energy: ##EQU8## where M  m2SH  ; weight of the secondary superheater metal. 
     C  m2SH  ; specific heat of the secondary superheater metal. 
     θ g2SH  ; temperature of the external gas of the secondary superheater. 
     α gm2SH  ; rate of heat transfer from the external gas of the secondary superheater to the metal 
     
         =α.sub.gm2SH, r.sup.(F.sub.gBF /F.sub.gBF, r).sup.0.8 (26) 
    
     α gm2SH , r ; rate of heat tranfer from the gas to the metal in a rated state. 
     F gBF  ; flow rate of the boiler gas. 
     F gBF , r ; rated gas flow rate of the boiler. 
     In addition, the gas temperature θ g2SH  of the secondary superheater becomes: ##EQU9## where H u  ; calorific power of the fual. 
     F f  ; flow rate of the fuel. 
     H a  ; ethalpy of the air. 
     F a  ; flow rate of the air. 
     H grf  ; enthalpy of the recirulation gas. 
     F grf  ; flow rate of the recirculation gas. 
     C Pg  ; specific heat of the gas. 
     K 1 , K 2  ; constants. 
     Here, when Equation (23) is arranged by putting the specific heat at constant pressure to be C p  =[∂H/∂θ] p , the following equation is obtained: 
     
         dx.sub.1 /dt=A.sub.11 x.sub.1 +A.sub.12 x.sub.2 +B.sub.11 u.sub.1 (28) 
    
     where 
     x 1  =θ s2SH   
     X 2  =θ m2SH   
     U 1  =θ s1SH   
     A 11  =-[C p  F s  +A s2SH  α ms2SH , r (F s2SH  /F s2SH ,r) 0 .8 ]/(V s2SH  ·γ s2SH  ·C p ) 
     A 12  : A s2SH  ·α ms2SH , r (F s2SH  /F s2SH , r) 0 .8 /(V s2SH  ·γ s2SH  ·C p ) 
     B 11  : C p  F s2SH  /(V s2SH  ·γ s2SH  ·C p ) 
     When Equation (25) is similarly generalized, the following equation is obtained: 
     
         dx.sub.2 /dt=A.sub.21 x.sub.1 +A.sub.22 x.sub.2 +B.sub.22 u.sub.2 (29) 
    
     where 
     u 2  =θ g2SH   
     A 21  =A s2SH  ·α ms2SH , r (F s2SH  /F s2SH , r) 0 .8 /(M m2SH  C m2SH ) 
     A 22  =-[A s2SH  ·α ms2SH , r (F gBF  /F gBF , r) 0 .6 +A s2SH  ·α m2SH , r (F s2SH  /F s2SH , r) 0 .8 ]/(M m2SH  ·C m2SH ) 
     B 22  =A s2SH  ·α gm2SH , r(F gBF/  F gBF , r) 0 .6 /(M m2SH  ·C m2SH ) 
     Equations (28) and (29) are state equations indicative of the characteristics of the secondary superheater. As the secondary-superheater gas temperature θ g2SH , the aforecited empirical equation (27) is used. 
     When the state transition equation of the secondary superheater is given by the following as described before, 
      (i+1)=φ(i) (i)+H(i) (i)                               (30) 
     the state transition matrix φ(t) and the drive matrix H(t) are given by the following equations: ##EQU10## Substituting Equations (28) and (29) into Equation (30), ##EQU11## where ##EQU12## By introducing them into Equations (5)-(8), the [Karman]--Kalman--filter can be constructed. 
     In this manner, the boiler system is put into the numerical expression model by the Kalman filter, whereby the most probable presumed value and predictive value of the main steam temperature can be obtained. Since, in the embodiment of this invention, the main steam temperature can be precisely measured actually, the heat transfer rates α mS2SH  and α gm2SH  included in the state transition matrix φ(i) are corrected and also the calorific power H u  included in the term of the drive matrix H(i) is corrected on the basis of the actually-measured values and predictive values and by the method explained with reference to FIG. 3, thereby permitting the presumption and prediction at higher precision. Of course, since the heat transfer rates are also included in the drive matrix H(i) as apparent from Equations (32) and (33), these parameters in H(i) are corrected in the aforecited corrections. 
     Now, the mill system will be described. In the mill system, it is essentially necessary to correctly presume the quantity of coal to be supplied to the burner. Since, however, the quantity of coal feed cannot be precisely measured, the quantity of primary air and the mill differential pressure (the difference between the pressure of the feed air to the mill and that of the air inside the mill) which can be measured comparatively precisely are also presumed and are used for model corrections so as to obtain a desired coal feed presumed value. 
     FIG. 5 is a constructional view of a simplified mill model. Here, the following equation holds: 
     
         A(dh.sub.c /dt)=F.sub.cf +F.sub.cr -F.sub.cg               (35) 
    
     where 
     F cg  =β g  N m   2  h c   
     F cr  =α g  F cg  =α g  β g  N m   2  h c   
     
         ∴A(dh.sub.c /dt)=F.sub.cf +(αβN.sub.m.sup.2 -βN.sub.m.sup.2)h.sub.c                              (36) 
    
     F CF  : quantity of coal feed to the mill. 
     h c  : height of coal in the mill. 
     F cg  : quantity of coal to be bitten between the balls and the table. 
     F cr  : quantity of coal to be recirculated. 
     β g  : rate of biting coal between the balls and the table. 
     α g  : rate of recirculation. 
     N m  : table revolution speed (constant). 
     A: table area. 
     In addition, 
     
         V=dγ.sub.c /dt=F.sub.cm -F.sub.cb                    (37) 
    
     where 
     F cm  =(1-α g )F cg  =(1-α g )β g  N m   2  h c   
     F cb  =γ c  F a   
     ∴V(dγ c  /dt)=(1-α g )β g  N m   2  h c  -F a , r γ c  -γ c , r F a   
     Kγ c  F a  =ΔP m   
     dΔP m  /dt=KF a , r dγ c  /dt (where K: constant between i and (i+1)) ##EQU13## where F CM  : quantity of coal staying in the mill. 
     F cb  : quantity of coal to be supplied to the burner. 
     γ c  : density of coal in the mill. 
     F a  : quantity of air to be supplied to the mill. 
     ΔP m  : mill differential pressure (air pressure at the mill inlet--air pressure in the mill). 
     γ c , r : γ c  during a rated operation. 
     F a , r : F a  during the rated operation. 
     V: internal capacity of the mill. 
     By arranging Equations (36) and (38), ##EQU14## It will be understood that, by bringing Equations (39) and (40) into correspondence with Equations (29) and (28), the Kalman filters can be correspondingly disposed also in the mill system likewise to the foregoing. For clarity, the corresponding relations of the variables in both the systems are listed below. 
     x 1  : ΔP m  - - - θ S2SH   
     x 2  : h c  - - - θ m2SH   
     u 2  : F cf  - - - θ g2SH   
     u 1  : F a  - - - θ S1SH   
     In an embodiment of this invention, the correction values of the state transition matrix φ(i) and the drive matrix H(i) as obtained by the foregoing method and noticed and are used for correcting also the gain of the feedback control system of the plant, whereby the control system is improved. It is now supposed that the feedback control system is a velocity type proportional-plus-integral control system. A computation equation is as given below: 
     
         P(i)=K.sub.I (i)·e(i, n)+K.sub.p (i){e(i, n)-e(i-1, n)}(38) 
    
     where 
     P(i): proportional-plus-integral control output at the i-th sampling. 
     K I  (i): integral control gain at the i-th sampling. 
     K p  (i): proportional control gain at the i-th sampling. 
     e(i, n): predictive derivation of the main steam temperture after n samplings at the time of the i-th sampling=T sp  (i, n)-T sms , p (i, n) 
     T sp  (i, n): desired value of the main steam temperature after n samplings at the time of the i-th sampling. 
     T sms , p (i, n): predictive value of the main steam temperature after n samplings at the time of the i-th 
     Assuming that the variation ΔT sms  of the main steam temperature versus the variation ΔF t  of the fuel flow rate conforms with the first-order lag characteristic (Equation (42)), the optimum control gains K p  (i) and K I  (i) are given as follows on the basis of the gain K g  (i) and time constant T g  (i) of the process at this time: ##EQU15## where s: Laplace operator. 
     β: constant. 
     K g  (i)=B 22  A 12  (H u  -K 1  K 2  F f   k2-1 )/[C PG  F GBF  (A 11  A 22  -A 12  A 21 )] 
     T gv  (i)=-(A 11  +A 22 )/(A 11  A 22  -A 12  A 21 ) 
     Accordingly, K g  (i) and T g  (i) may be corrected similarly to the corrections of φ(i) and H(i). 
     To sum up, FIG. 6 illustrates the concept of the control blocks in the case of the application of this invention to the thermal power plant. Numeral 601 designates a portion for generating an output desired-value signal, which is converted into a coal feed quantity request signal by a function generator 603. Shown at numeral 605 is an adder, which adds the output of the function generator 603 and a correction signal given from an adder 627. Numeral 607 indicates an adder, which adds the output of the adder 605 and a coal feed quantity presumed value supplied from a mill model 615 through a line 616. Numerals 609 and 611 indicate proportional-plus-integral control portions, which are respectively supplied with correction signals from the model 615 through lines 610 and 612 so as to establish the gains given by Equation (43). An adder 613 adds the outputs of the control portions 609 and 611. Numeral 30 indicates the mill system, and the quantity of coal feed to this mill system is controlled by the output of the adder 613. The various data of the mill system are sampled into the mill model 615 through lines 351 and 352. Fuel is supplied from the mill system 30 to the boiler system 40. The various data of the boiler system are sampled into a boiler model 617 through a line 451. The coal feed quantity presumed value of the mill model is also sampled into the boiler model through the line 616. The main steam temperature predictive value after a predetermined time as obtained from the boiler model is compared by an adder 621 with a desired value after the predetermined time as given by a main steam temperature desired-value signal generator portion 619, and the difference is applied to proportional-plus-integral control portions 623 and 625. These control portions are supplied with correction signals from the model 617 through lines 624 and 626 so as to establish the gains given by the foregoing equation (43). The outputs of the respective control portions 623 and 625 are applied to the adder 627, the output of which is added to the coal feed quantity request signal as stated above and corrects this signal so as to eliminate an error in the main steam temperature. 
     FIG. 7 is a waveform diagram for explaining the effect of this invention, and exemplifies a heating control in a thermal power plant. In the figure, numeral 71 designates a desired temperature rise curve, numeral 72 the curve of a control result according to the embodiment of this invention, and numeral 73 the curve of a control result according to an example of the U.S. patent cited before. As apparent from the comparison between the curves 72 and 73, the embodiment of this invention produces a better control result.