CN1570035A - Catalytic reforming 17 lump reaction model modelling approach - Google Patents

Abstract

The invention discloses a catalytic reforming 17 lump reaction model modelling approach characterized in that, (1) the reaction feedstock is divided into 17 lumped components, (2) recapitalization reaction network containing 17 reactions is included, (3) method for determining model parameter is provided, (4) continuance reforming catalyst coke and deactivation model is created, (5) differentiated accurate one-dimensional searching algorithm is introduced into the classical BFGS alogorithm, which is a parameter estimation.

Description

Catalytic reforming 17 lumpedreaction model modeling method

Technical Field

The invention relates to a modeling method of a catalytic reforming 17 lumped reaction model.

Background

The catalytic reforming device is one of important secondary processing devices in the petroleum processing process and one of complex catalytic reaction systems. The chemical reaction generated in the reforming process mainly comprises dehydrogenation and isomerization of naphthene,Five types of reactions such as alkane dehydrocyclization, arene hydrogenolysis, alkane hydrocracking and the like. The reforming feed is typically refined naphtha (primarily C)₆～C₁₁Alkanes, cycloalkanes and aromatics) with nearly 300 pure components, the reaction being more countless. The reforming device not only can produce high-octane gasoline, but also can provide a large amount of aromatic hydrocarbons as important chemical raw materials, is an indispensable important production device for modern oil refineries and petrochemical plants, and is one of important sources of economic benefits of petrochemical enterprises.

Because the whole production equipment is old, the process technology and the automation level are laggard, and the like, the reforming device in China cannot compete with foreign countries in the aspects of energy consumption, average production cost and oil quality. This is very disadvantageous for the intense international market competition directly after the addition of WTO, and therefore how to improve the device economy has become one of the important concerns of the enterprise decision-making layer. The method for guiding production design, optimizing device operating conditions, eliminating production bottleneck, realizing process monitoring and fault diagnosis by utilizing the reaction mechanism model of the reforming device is one of the most effective methods for improving the economic benefit of the device.

Due to the importance of the reforming model, a great deal of research and development work is carried out on the reforming model by many researchers in domestic and foreign research institutions and universities. As described above, the reforming reaction involves many components and reactions, and therefore, the research is generally conducted by a lumped method. Smith [1]1959]The four lumped reforming reaction model was first proposed, i.e. the naphtha feed was divided into alkanes, cycloalkanes, aromatics and cracked products (C)₅-) and the like are lumped (contains no hydrogen component, the same applies below) and only 4 main chemical reactions of alkane cyclization, cycloalkane dehydrogenation, alkane cracking, and cycloalkane cracking are considered. At the end of the 70 s, Zingiberam officinale (2) of the institute of great Liang Hua in China]The catalytic reforming reaction system is researched, hydrocarbons with different carbon atoms are distinguished, and a 14-lumped model containing 26 reactions is provided, so that the catalytic reforming reaction system is greatly improved compared with a Smith model. Ramage et al, Mobil, 80 s [3]]A13 lumped kinetic model is proposed comprising 24 reactions, which also distinguishes between hydrocarbons of different carbon numbersAnd pentabasic, hexahydric, naphthenic hydrocarbons, without distinguishing the n-isoparaffins. At the end of the 80 s, Fromert et al [4]]Not only considers the different properties among normal isoparaffin and five-membered six-membered cycloalkane, but also considers C₉The hydrocarbons have a comparable composition in the reforming feed and the cracked product C₅Finer resolution, a 28 lumped model is proposed containing 84 reactions. In the early 90 s of the year,weng Hui Xin, Sun Shao Zhuang and other people of China east university of science [5]On the basis of 13 lumped models proposed by Ramage, the normal isoparaffin is further distinguished, and a 16 lumped reforming reaction model containing 27 reactions is developed, but the model does not distinguish C₈The above hydrocarbons. Meanwhile, Luoyang petrochemical engineering Co., Ltd. [6]deazan et al]A reforming kinetics model similar to the 28-lumped model of fromet was also developed, but varied in the specific partitioning of components and number of reactions, which contained 68 reactions. Other similar models have been developed in recent years, such as Jorge [7]]24 lumped, contains a reforming kinetics model of 71 reactions, etc.

However, few detailed reports or publications of currently developed reforming models have been made or reported for reasons such as commercial security.

In general, the disadvantages of the existing reforming reaction models are:

1. lumped partitioning does not meet model accuracy or industrial analysis requirements. Like the 4-lumped model of Smith, the lumped division is too little, the difference of the reaction performance among hydrocarbon components with different carbon numbers is ignored, and the model precision is not high. Other reaction models distinguish not only hydrocarbon components with different carbon numbers, but also normal isoparaffin or quinary-senary cycloalkane, and the lumped division is too complex. At present, most of the reformers in China have limited laboratory analysis conditions, and generally only C in reformed raw oil or product oil is analyzed₁～C₉₊The composition of the families of alkane, cycloalkane and arene can not meet the analysis requirements of the model, and brings difficulty to the industrial application of the model.

2. The model is too complex and the parameter estimation is difficult. Particularly, most of models developed since the 90 s have more than 20 lumped numbers, even 64 lumped numbers, and the generated results are that the number of reactions is increased sharply, the dimensions of the model equation and the mechanism parameter variable to be estimated are too high, and the parameter estimation algorithm is difficult to converge.

Aiming at the two main defects, the invention provides a simplified mechanism model which has more reasonable lumped division and less reaction number and covers the main reforming reaction based on an actual industrial reforming device and fully considers the limited laboratory analysis conditions of most of the reforming devices in China at present. The model distinguishes C₆～C₉₊And the group composition of alkanes, cycloalkanes and aromatics, and the refining of the cracked product to C₁～C₅There were 17 lumped components (excluding the hydrogen component) comprising only 17 reactions. The model greatly reduces the requirements on raw material analysis, solves the technical problem that the parameter estimation algorithm is difficult to converge due to excessive component aggregation and high dimension of the parameter to be estimated, and greatly facilitates the practical industrial application of the model.

Reference documents:

[1]Smith，R.B.(1959)Kinetic Analysis of Naphtha Reforming with PlatinumCatalyst，Chem.Eng.Prog.，55(6)：76～80.

[2]study of catalytic reforming reaction kinetics [ D]zonula (1978): large linkage of the institute of chemico-physical research.

[3]Ramage，M.P.，et al.(1980).Dev.of Mobil kinetic reforming model.Chem.Eng.Sci.，35(1)：41.

[4]Froment，G.F.The Kinetic of Complex Catalytic Reactions，Chem.Eng.Sci.，1987，42(5)：1073～1087.

[5]Kung hui, sun shazhuang et al (1994). establishment of catalytic reforming lumped dynamics model (I), chemical bulletin, 45 (4): 407 to 411.

[6]Dean et al (1995) establishment of catalytic reforming reaction kinetic model and its industrial application (1) refinery design, 25 (6): 49 to 51.

[7]Jorge A.J.，Eduardo V.M.(2000).Kinetic Modeling of Naphtha CatalyticReforming Reactions.Energy&Fuels，14(5)：1032～1037.

Disclosure of Invention

The invention aims to provide a modeling method of a catalytic reforming 17 lumped reaction model. The modeling method is characterized in that:

(1) subdivision of the reaction feed into 17 lumped components

a) Alkane: methane (C)₁) Ethane (C)₂) Propane (C)₃) N-isomeric butanes (C)₄) N-isomeric Pentanes (C)₅) N-isomeric Hexane (P)₆) N-isomeric Heptane (P)₇) N-iso-octane (P)₈) And an alkane (P) having 9 or more carbon atoms₉₊) Equal to 9 lumped;

b) cycloalkane: methylcyclopentane and cyclohexane (N)₆) Cycloalkane having 7 carbon atoms (N)₇) Cycloalkane having 8 carbon atoms (N)₈) Cycloalkane (N) having 9 or more carbon atoms₉₊) Equal 4 lumped;

c) aromatic hydrocarbons: benzene (A)₆) Toluene (A)₇) Xylene and ethylbenzene (A)₈) And heavy aromatic hydrocarbons containing 9 or more carbon atoms (A)₉₊) Equal 4 lumped;

(2) reforming reaction network comprising 17 reactions

The reaction network comprises 4 alkanes (P)₆，P₇，P₈，P₉₊) Dehydrocyclization to give cycloalkane reaction, 4 cycloalkanes (N)₆，N₇，N₈，N₉₊) Dehydrocyclization to generate aromatic hydrocarbon reaction, toluene hydrogenolysis to generate benzene and methane, xylene hydrogenolysis to generate toluene and methane, trimethylbenzene hydrogenolysis to generate xylene and methane, and methylethylbenzene hydrogenolysis to generate toluene and ethane, 4 aromatic hydrocarbon hydrogenolysis reactions, 5 alkane hydrocracking reactions and other 17 reforming reactions;

(3) method for determining model parameters to be estimated

Estimating only 17 corrected frequency factor parameters in the mechanism model;

(4) continuous reforming catalyst coking and deactivation model

The integral term in the catalyst coking model adopts a sectional integral method, and the inactivation model introduces a reactor average activity function, namely the catalyst activity of a single reactor is considered to be the same;

(5) parameter estimation algorithm

A differential accurate one-dimensional search algorithm is introduced on the basis of a classic BFGS algorithm, namely a Hessian matrix of a target function is approximated by utilizing the product of a Jacobian matrix of the target function, so that an analytic expression of the optimal step length of one-dimensional search is solved.

The invention has the advantages that:

1) the lumped division not only basically accords with the dynamic mechanism of the reforming reaction, but also accords with the analysis conditions of industrial reforming devices in China, and the industrial application is convenient;

2) the model parameters needing to be estimated are few, and the estimation difficulty is low;

3) catalyst coking and deactivation models are integrated;

4) the parameter estimation algorithm has superior convergence performance.

Drawings

FIG. 1 is a schematic flow diagram of a catalytic reformer;

FIG. 2 is a schematic diagram of a 17 lumped reforming reaction network;

FIG. 3 is a schematic diagram of a reforming radial reactor configuration in which: the device comprises a catalyst inlet 1, an oil gas inlet 2, a catalyst bed layer 3, a reactor cylinder 4, a catalyst outlet 5 and an oil gas outlet 6;

FIG. 4 is a schematic diagram of a catalyst bed structure.

Detailed Description

The terms of the present invention are explained:

catalytic reforming: is a secondary processing process of petroleum. I.e. naphtha (containing C)₆～C₁₁Hydrocarbon substance) in the presence of a catalyst, and performing chemical reactions such as aromatization, isomerization and the like of molecular structure rearrangement to prepare aromatic hydrocarbon or high-octane gasoline.

Lumped: for reactants containing a plurality of pure components, a compound with similar kinetic properties and physicochemical properties can be regarded as a virtual component, so as to achieve the purpose of reducing the number of the reactants. This method of component clustering we call aggregation.

Aromatic hydrocarbon: the quality index for evaluating the quality of the reforming raw material comprises the following calculation methods:

$A_r\%=0.93C_6^N\%+0.94C_7^N\%+0.95C_8^N\%+\mathrm{\Σ}{A}_r^0\%$

wherein C is₆ ^N％，C₇ ^N％，C₈ ^N% represents the content of cycloalkanes containing 6, 7 and 8 carbon atoms respectively in percentage by mass;

representing the total aromatics content of the feed.

Hourly space velocity of reaction liquid: is the feed flow (m) of the reforming liquid³H) and reforming catalyst loading (m)³) The ratio of (a) to (b) indicates the operating load of the reforming reactor.

1. Description of catalytic reforming apparatus

The catalytic reforming device is one of important secondary processing devices in the petroleum processing process, not only can produce high-octane gasoline, but also can provide a large amount of aromatic hydrocarbons which are used as important chemical raw materials. The reaction separation part of the continuous reforming device (shown in figure 1) mainly comprises 4 overlapped reforming reactors, 1 catalyst regenerator, 1 four-in-one heating furnace, 1 high-pressure separation tank and 1 hydrogen purification system. Wherein the reforming reactor is a core component of the apparatus. The process flow of the catalytic reformer shown in fig. 1 is: the reforming feed is firstly mixed with circulating hydrogen from a high-pressure separation tank, exchanges heat with a reforming reaction product (reforming generated oil), enters a heating furnace to be heated to a certain reaction temperature and then enters a first reforming reactor to carry out reforming reaction, and generated oil gas enters the heating furnace to be heated to a certain reaction temperature after coming out of the first reforming reactor and then enters a second reforming reactor to carry out reforming reaction until reaching a fourth reforming reactor. And the high-temperature reformed oil gas from the fourth reforming reactor exchanges heat with the reformed feed, enters an air cooler for further cooling, and then enters a high-pressure separation tank for gas-liquid separation. And one part of the separated hydrogen-containing gas is used as circulating hydrogen, and the other part of the separated hydrogen-containing gas enters a hydrogen purification system for purification and then is sent out of the device to be used as a reforming byproduct. And during the reforming reaction, the reforming catalyst slowly moves from top to bottom in the 4 reactors, the coked catalyst is sent into a catalyst regenerator for regeneration after coming out of the fourth reforming reactor, and the regenerated catalyst is sent into the first reforming reactor as a fresh catalyst to continue the reforming reaction, so that the process is circularly carried out.

(2) Reforming reaction model establishment

a) Lumped partitioning and reaction network

The catalytic reforming process is one of complex catalytic reaction systems, and the generated chemical reactions mainly comprise five types of reactions such as naphthene dehydrogenation, isomerization, paraffin dehydrocyclization, arene hydrogenolysis, paraffin hydrocracking and the like. The reforming feed is typically refined naphtha (primarily C)₆～C₉₊Alkanes, cycloalkanes and aromatics) with nearly 300 pure components, the reactions of which are more numerous and must be modeled using a lumped method. Considering that many domestic reformers cannot analyze very fine compositions (typically only C is analyzed)₁～C₉₊Group composition of alkane, cycloalkane and arene) or the analysis task is too heavy and the analysis precision is not high, so that a large amount of industrial data samples meeting the analysis requirements cannot be provided for model parameter estimationCounting; on the other hand, too much aggregation causes too high dimension of the model equation and the parameter to be estimated, and also causes the problem that the parameter estimation algorithm is difficult to converge. Therefore, although the more the number of component lumped, the higher the model accuracy is likely to be, the more difficult it is to be industrially applied. In consideration of these practical situations and the reforming reaction mechanism, we do not distinguish the difference of the reforming reaction performance of the normal isoparaffin and the five-membered six-membered cycloalkane, but rather, the C₆～C₉₊Respectively, and subdividing the cracked gas into C₁～C₅Alkanes, and thus the gas composition in the recycle hydrogen and the by-product hydrogen can be calculated. The detailed lumped division is as follows:alkane: methane (C)₁) Ethane (C)₂) Propane (C)₃) N-isomeric butanes (C)₄) N-isomeric Pentanes (C)₅) N-isomeric Hexane (P)₆) N-isomeric Heptane (P)₇) N-iso-octane (P)₈) And an alkane (P) having 9 or more carbon atoms₉₊) Equal to 9 lumped;

cycloalkane: methylcyclopentane and cyclohexane (N)₆) Cycloalkane having 7 carbon atoms (N)₇) Cycloalkane having 8 carbon atoms (N)₈) Cycloalkane (N) having 9 or more carbon atoms₉₊) Equal 4 lumped;aromatic hydrocarbons: benzene (A)₆) Toluene (A)₇) Xylene and ethylbenzene (A)₈) And heavy aromatic hydrocarbons containing 9 or more carbon atoms (A)₉₊) Etc. 4 lumped.

We thus obtained a simplified reforming reaction model with 17 lumped components (excluding the hydrogen component) and containing only 17 reactions. The lumped reactive network is shown in figure 2.

The reaction network in fig. 2 includes 17 reforming reactions such as 4 alkane dehydrocyclization to generate cycloalkane reaction, 4 cycloalkane dehydrocyclization to generate aromatic hydrocarbon reaction, 4 aromatic hydrocarbon hydrogenolysis reaction, and 5 alkane hydrocracking reaction, which are specifically as follows:

alkane dehydrogenation reaction (reversible):

P_i_N_i+H₂，i＝6，7，8，9+

cycloalkane dehydrogenation reaction (reversible):

N_i_A_i+3H₂，i＝6，7，8，9+

alkane hydrocracking reaction:

$C_5+H_2\→1/2\·\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{C}_i$

$P_5+H_2\→1/3\·(\underset{i=1}{\overset{5}{\mathrm{\Σ}}}{C}_i+C_3)$

$P_5+H_2\→1/3(P_6+\underset{i=1}{\overset{5}{\mathrm{\Σ}}}{C}_i)$

$P_8+H_2\→1/4(P_7+P_6+\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{C}_i+C_4)$

$P_{9+}+H_2\→1/4(P_8+P_7+P_6+\underset{i=1}{\overset{5}{\mathrm{\Σ}}}{C}_i)$

and (3) carrying out hydrogenolysis reaction of aromatic hydrocarbon:

，

，

in building this reaction network, we consider the following in terms of aggregation and reaction mechanism:

n-isoparaffins are not distinguished. Previous studies of reforming reactions have shown that the differences in the reaction performance of n-isoparaffins are essentially negligible, so we have not differentiated n-isoparaffins in this model in order to alleviate the task of industrial analysis. Research has shown that, for reformate octane number estimates, despite differences in the octane number of n-isoparaffins, for the same reformate feed, changes in octane number due to changes in isomerization reaction resulting from changing operating conditions are essentially negligible and the reformate octane number is approximated as a quasi-linear function of the aromatic content [1]]。

Five-membered, six-membered cycloalkane is not distinguished.Since the conversion rates of five-membered hexatomic cycloalkanes (except methylcyclopentane and cyclohexane) in the reforming reaction are all very similar, and the isomerization reaction is fast, we do not distinguish five-membered hexatomic cycloalkanes.

Distinction C₈And C₉₊And (4) group composition. Due to C₉₊The group composition is a significant proportion of the reforming feed (about 30% mass fraction), and considering that the reaction model must yield the triphenyl (benzene, toluene, xylene) when used in an aromatics production scheme, C must be substituted₈From C₈₊Is separated out.

The cracked product is subdivided. To obtain a more accurate recycle hydrogen purity, and specific composition of Fuel Gas (FG), Liquefied Petroleum Gas (LPG) and stabilized gasoline, we subdivide the cracked products into C₁～C₅An alkane.The benzene ring only carries out dealkylation reaction (hydrogenolysis reaction) and does not carry out hydrocracking reaction on the benzene ring, and the generated micromolecular alkane is mainly C₁～C₂It is considered that ethylbenzene and propylbenzene generated in the reaction are negligible, and the probability of hydrogenolysis of each alkyl group having a carbon number is uniform without simultaneously removing a plurality of alkyl groups.

Since the hydrocracking process of the cycloalkane is ring-opening to generate the alkane, we mainly consider the hydrocracking of the alkane, and consider that all the hydrogenation reactions occur with equal probability to simplify the reaction model.

b) Catalyst coking and deactivation model

The reforming catalyst can generate coking phenomenon on the surface during the operation process, and the deposited coke can block the pores of the catalyst and cover the metal/acid active center of the catalyst, thereby causing the catalyst to be deactivated. In the reforming reaction, the reactants, products and many intermediate molecules are likely to become coke precursors or combine with each other and condense into high molecular weight deposits, so that the catalyst is prone to carbon deposition and the mechanism is complicated. In order to investigate the reforming reaction characteristics under the catalyst deactivation condition and accurately predict the reformate yield and coke content in the operation process of an industrial plant to optimize the reforming reaction conditions and the catalyst regeneration conditions, we propose the following catalyst coking and deactivation models:

φ＝exp(-α·CK(t)) (2)

in the above formula r_c(t) is the catalyst coking rate, kg/kg cat/s; r is_cThe catalyst has the advantages of high activity, high catalyst activity and high catalyst coking rate.

For convenience of model solution, for r_cWe propose the following correlation

Wherein k is a constant, and a, b, c, d and e are variable coefficients; p_hIs hydrogen partial pressure (Mpa), A_rWt% for feeding aromatic hydrocarbon; t is_FEPFeed end point, deg.C; t is₀The temperature is 165 ℃; t is the reaction temperature, DEG C. Obviously, the above formula is easily converted into a linear regression model through variable replacement, so that the coefficient can be easily regressed. For simplificationSee, we get r_cThe temperatures and pressures in DEG are the reactor mean inlet temperature (WAIT, DEG C) and the mean pressure (1 arithmetic mean of the pressure at the inlet and the pressure in the upper partial tank, MPa), respectively, so that r for all reforming reactors_cEqual and are denoted as r_cDegree. Substituting the formula (2) into the formula (1) to obtain a model equation of the catalyst coking content (kg/kg. cat) which varies with the running time t of the device

Wherein the catalyst activity function φ is calculated from the formula (2).

For the continuous reforming device, the catalyst is continuously regenerated, namely, the regenerated fresh catalyst enters from the first reactor and slowly passes through each reactor under the action of gravity, reactant gas passes through gaps of catalyst particles to generate reforming reaction in contact, and the spent catalyst continuously flows out from an outlet of the fourth reactor and enters into the regenerator for regeneration. It is clear that the carbon content of the catalyst is a function of the residence time of the catalyst in the reactor or the distribution of the height of the catalyst bed. Considering the reactor characteristics of a continuous reformer, we can obtain the following dynamic model equation of the coke content (kg/kg cat) of the catalyst as a function of the operating time t and the height z of the catalytic bed,

$t_z=t_{1j-1}+\frac{{\mathrm{\Ω}}_j\·{\mathrm{\ρ}}_c}{v_c}(z-H_{1j-1}),t_{10}=0,t_{1j}=\underset{k=1}{\overset{j}{\mathrm{\Σ}}}\frac{{\mathrm{\Ω}}_k\·{\mathrm{\ρ}}_c}{v_c}{H}_k$

$H_{10}=0,H_{1j}=\underset{k=1}{\overset{j}{\mathrm{\Σ}}}{H}_k,j=1,...,n$

and catalyst coking rate (kg/h) equation

Coke＝0.95×CK/(1-0.95×CK)·v_c(6)

Middle omega of the above formula_jThe cross section of the catalytic bed layer of the jth reactor, m²；ρ_cAs catalyst density, kg/m³；z，H_k，H_1jOn an axial scale, the height of the catalytic bed of the kth reactor, from the 1 st reactor to the jth reactorThe total height of the catalytic bed layer of the reactor, m; v. of_cKg/h is the catalyst circulation rate. Generally, the carbon content in the catalyst coke is about 95%.

To the coke content at the outlet of the last reactor

In the above formula r_cDegree (t) is the curve of coking rate with fresh catalyst as a function of unit run time t, t_cThe catalyst cycle regeneration time. Obviously, only r needs to be known_cThe coking content can be conveniently estimated by integrating the value of the curve by the formula. Since the operating condition fluctuation is not so large, if the integration time is divided in hourly, r_cThe degree (t) being substantially constant during each integration period, whereby fractional numerical integration is usedThe accuracy will be relatively high, i.e. the integral term can be expressed as

Substituting the formula (5) into the formula (2) to obtain the catalyst activity function

For simplicity, we assume that the catalyst moves uniformly within the reactor and that the catalyst activity is the same for the same height of the catalyst bed, ignoring radial variations in the coke content of the catalyst in the reactor. As can be seen from the radial reactor configuration (fig. 3, 4), reformate almost simultaneously traverses the longitudinal cross-section of the reactor, which means that catalysts of different activity participate in the reaction for the same residence time (with the same radial position). If we divide the reaction oil gas into N equal parts along the longitudinal bed, the material balance equation of the reforming model can be rewritten as

Arranged in the above formula to obtain

$\frac{d\stackrel{\→}{Y}}{\mathrm{dR}}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\φ}}_i/N\·\frac{2\mathrm{\πR}\·H}{\mathrm{LHSV}\·V_c}{K}_r\·\stackrel{\→}{Y}=\stackrel{\‾}{\mathrm{\φ}}\·\frac{2\mathrm{\πR}\·H}{\mathrm{LHSV}\·V_c}{K}_r\·\stackrel{\→}{Y}----\left(11\right)$

Where phi is a function of the average activity,

$\stackrel{\‾}{\mathrm{\φ}}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\φ}}_i/N,{\mathrm{\φ}}_i={\mathrm{\φ}}_i(z,t).$

we assume that φ (z, t) is a continuous function, then according to the principle of calculus, it can be described by Riemann integration when N → ∞ $\stackrel{\‾}{\mathrm{\φ}}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\φ}}_i/N,$ Namely have

$\stackrel{\‾}{\mathrm{\φ}}\left(t\right)=\underset{N\→\∞}{\lim }\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\φ}}_i/N=\∫\mathrm{\φ}(z,t)\mathrm{dz}/H----\left(12\right)$

For any reactor, the catalyst average activity function can be expressed as

${\stackrel{\‾}{\mathrm{\φ}}}_j\left(t\right)={\∫}_{h_j}^{h_j+H_j}\mathrm{\φ}(z,t)\mathrm{dz}/H_j,j=1,...,n----\left(13\right)$

Wherein h is_j＝H_1j-1，H_jRespectively the initial value of the axial coordinate of the jth reactor (scaled down from the top) and the overall height, m. The average catalyst activity in each reactor at time t is expressed by substituting equation (9) into the above equation (integration time interval is set to 1h)

It follows that the average activity function of the catalyst in the same reactor can be considered the same and is only a function of the reactor and the run time t. For different reforming reactions, the catalyst activity function is considered to be the same, i.e. to have

c) Equation of model

Strictly speaking, the reforming reactions are complex heterogeneous catalytic reactions. To simplify the complexity of the model, we have used the results of previous studies that suggest that each reforming reaction is a simple first order reaction with respect to the hydrocarbon composition and is exponential in hydrogen partial pressure. Because the hydrogen-hydrocarbon ratio is large and the hydrogen partial pressure in the reactor does not change greatly, the whole reaction system can be regarded as a pseudo-first-order homogeneous reaction system. Based on this theory, and in conjunction with the previous reaction network, we derive 17 reaction rate equations as follows:

alkane dehydrogenation reaction (reversible): r is_j＝dY_Pi/d(1/LHSV)＝k_j·(Y_Pi-Y_Ni/K_epj)，j＝1～4

Cycloalkane dehydrogenation reaction (reversible): r is_j＝dY_Ni/d(1/LHSV)＝k_j·(Y_Ni-Y_Ai/K_epj)，j＝5～8

Alkane hydrocracking reaction: r is_j＝dY_Pi/d(1/LHSV)＝k_j·Y_Pi，i＝5～9+，j＝9～13

And (3) carrying out hydrogenolysis reaction of aromatic hydrocarbon: r is_j＝dY_Ai/d(1/LHSV)＝k_j·Y_Ai，i＝7～9+，j＝14～17

Wherein Y is the molar flow of each lumped component, kmol/h; LHSV is liquid hourly space velocity, 1/h; k_epjIs the equilibrium constant of the reversible reaction; k is a radical of_jIs a reaction rate constant which is a function of the reaction temperature T, the reaction pressure P and the hydrogen-oil molar ratio r and is expressed by

$k_j=k_{0j}\·\exp (-E_j/\mathrm{RT})\·P_h^{{\mathrm{\θ}}_j}\·{\stackrel{\‾}{\mathrm{\φ}}}_i,i=\mathrm{1,2},...,n;j=\mathrm{1,2},...,17----\left(16\right)$

Wherein k is_0jAs a reaction frequency factor, h^-1·Mpa^-θj；E_j，θ_jRespectively representing reaction activation energy (KJ/mol) and pressure index; r is a gas constant of 8.314 J.mol^-1·K^-1；P_hIs the partial pressure of hydrogen in the reactants, P_h≈P·r/(r+1)，Mpa；φ_iIs a function of the average activity of the catalyst of the ith reactor and is calculated by a formula of a catalyst deactivation model (14); n is the number of reactors.

To simplify the model equation, we make the following assumptions based on the structural features of the radial reactor and the previous study [6]: (1) the axial temperature is uniform; (2) the axial concentration is uniform; (3) the reactor is considered an ideal adiabatic reactor. Based on the principle of material balance and energy balance, we can obtain the following reforming reaction model equation set (ordinary differential equation set with initial value),

whereinMolar flow of 18 lumped components (including hydrogen component), kmol/h; wherein R is the coordinate of the flow direction of the reaction oil gas, m; h is the reactor height, m; Δ H_jKJ/kmol as heat of reaction;the specific heat of 18 lumped components, KJ/kmol/K.

K_r＝C·K_vIn which K is_vIs a matrix of rate constants for each reaction, whose elements are k_iOr-k_i/K_epi，

Namely, it is

C is a matrix of stoichiometric coefficients in the reaction to which the components are involved, i.e.

(3) Solving of model equations

The catalytic reforming model equation is a typical initial value problem of a nonlinear ordinary differential equation set (ODE) according to the formula (17).

For general nonlinear ODE initial value problem

$d\stackrel{\→}{y}/\mathrm{dx}=\stackrel{\→}{f}(x,\stackrel{\→}{y}),\stackrel{\→}{y}\left(x_0\right)=\stackrel{\→}{y}----\left(18\right)$

Because of the complexity of the equations, it is difficult to solve the analytical solution, and a numerical solution method providing an initial value is generally used. At present, many numerical algorithms of ODE are mature, and are mainly divided into a single-step method and a multi-step method. The single-step method is the most common algorithm and is characterized by calculating y_i+1The formula containing only the previous step y_iTypical algorithms include euler's forecast-correction method, longgutta method), and the like. The multi-step method is also called as linear multi-step method and is characterized in that the left side and the right side of the structural formula are respectively linear expressions of y and fAnd calculate y_i+1The formula includes the first multiple steps y_i，y_i-1…, there are known algorithms such as Adams-Bashforth method (explicit), Adams-Moulton method (implicit), and Gear method. The advantage of the multi-step method is that the function is calculated

Has less times when

When the algorithm is complex, the calculation time can be obviously reduced, but the algorithm needs former multi-step information for starting. The high reliability, high precision and less calculation time of the algorithm are the criteria for evaluating the quality of the algorithm. In fact, since the reforming reaction rate is very different (for example, the dehydrogenation reaction rate of naphthene is almost more than 100 times of that of hydrocracking), the model equation has the problem of rigidity (stiffness), that is, the numerical algorithm may be unstable and the calculation result is very different by adopting different step sizes. It has been shown that explicit numerical algorithms are not likely a-stable (i.e., the algorithms will have stability regardless of step size), so explicit algorithms cannot be used for ODE solutions with rigid problems. Since the classical RK method is an explicit algorithm, it has irreconcilable contradictions in solving rigid ODEs: when the step length islarge, not only is the solution inaccurate, but also the algorithm is unstable; conversely, if the step size is small, rounding errors are introduced and computation time is greatly increased. According to the Taylor expansion idea, the invention provides an approximate analysis method with a simpler form, and the approximate analysis method is combined with an implicit Gear method, so that the problem of the rigidity of an equation is solved, the stability and the calculation precision of an algorithm are ensured, and the calculation time is shortened.

a) Approximate analytic numerical algorithm

For convenience, we change equation (17) to

Wherein $K_r^{\′}=\frac{2\mathrm{\π}\·H}{\mathrm{LHSV}\·V_c}{K}_r,K_r=K_r\left(T\right).$

From formula (19), if the step size Δ R is sufficiently small, the reaction temperature T is considered to be substantially constant, matrix K'_rApproximately as a constant matrix. For equation (19), itsgeneral mathematical expression can be written as

Wherein the content of the first and second substances,is an 18-dimensional function vector, and K is an 18 x 18 constant matrix. Assumption function

If there is a continuous derivative of any order in a domain at point x, then we can do this for a functionPerforming Taylor series expansion, i.e. by

$\stackrel{\→}{y}(x+\mathrm{\Δx})=\stackrel{\→}{y}\left(x\right)+\stackrel{\→}{y}{\left(x\right)}^{\′}\·\mathrm{\Δx}+\stackrel{\→}{y}{\left(x\right)}^n\·\mathrm{\Δ}{x}^2/2+\·\·\·=\underset{n=0}{\overset{\∞}{\mathrm{\Σ}}}\stackrel{\→}{y}{\left(x\right)}^{\left(n\right)}\·\mathrm{\Δ}{x}^n/n!----\left(21\right)$

From the formula (20), it is apparent that

$\stackrel{\→}{y}{\left(x\right)}^{\′}=x\·K\·\stackrel{\→}{y}\left(x\right),\stackrel{\→}{y}{\left(x\right)}^{\′\′}=K\·\stackrel{\→}{y}\left(x\right)+x^2\·K^2\·\stackrel{\→}{y}\left(x\right),\·\·\·$

The derivatives of each order are substituted into the formula (21) and are arranged to obtain

$\stackrel{\→}{y}(x+\mathrm{\Δx})=(I+K\·{\mathrm{\Δx}}^{\′}+K^2\·{\mathrm{\Δx}}^{\′2}/2+\·\·\·)\·\stackrel{\→}{y}\left(x\right)----\left(22\right)$

Where Δ x ═ (x + Δ x/2) · Δ x, and I is an 18 × 18 identity matrix. Accordingly, we can obtain an approximate analytical solution of equation (19) for material balance

$\stackrel{\→}{Y}(R+\mathrm{\ΔR})=(I+K_r^{\′}\·{\mathrm{\ΔR}}^{\′}+{K_r^{\′}}^2\·{\mathrm{\ΔR}}^{\′2}/2+\·\·\·)\·\stackrel{\→}{Y}\left(R\right)----\left(23\right)$

Wherein Δ R ═ (R +1/2 · Δ R) Δ R.

It can be shown that if the first 2 terms are taken on the right side of equation (22), the above equation becomes the well-known first-order euler method; an approximate improved Euler method taking the first 3 terms as second order; and the first 5 terms are the approximate RK method of fourth order (error term is ${\mathrm{\ΔR}}^{\′6}\·{(R+{\mathrm{\ΔR}}^{\′}/2)}^2\·K_r^{\′4}/24$ ) The more terms are taken, the higher the precision. In consideration of the influence of calculation time and accumulated errors, the first 5 terms of the equation (22), namely the fourth-order approximation RK method, are taken to solve the model equation. Obviously, equation (22) not only unifies the basic classical explicit algorithms, but also has a very simple form. Therefore, the approximate analytic numerical algorithm has obvious advantages in programming compared with the classical RK method.

b) Gear numerical algorithm

The Gear method is an implicit linear numerical algorithm, and the general expression of the algorithm is

$\underset{j=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\α}}_j\·{\stackrel{\→}{y}}_{n+j}=h\·{\mathrm{\β}}_k\·{\stackrel{\→}{f}}_{n+k}----\left(24\right)$

The above formula is also called k-step Gear method, wherein α_j，β_kAnd h is the step size. For the four-step Gear method, itThe algorithm equation is

${\stackrel{\→}{y}}_{n+4}-\frac{48}{25}{\stackrel{\→}{y}}_{n+3}+\frac{36}{25}{\stackrel{\→}{y}}_{n+2}-\frac{16}{25}{\stackrel{\→}{y}}_{n+1}+\frac{3}{25}{\stackrel{\→}{y}}_n=\frac{12}{25}h\·{\stackrel{\→}{f}}_{n+4}----\left(25\right)$

The solution of the initial value problem (20) using equation (25) requires 4 starting valuesExcept that

The remaining 3 initial step values are calculated by a one-step method except for the given step, and the four-order approximation analysis method is adopted for calculation in the invention. For implicit linear algorithms, this is generally due to

Is a about

Is thus solved from equation (25)Generally, an iterative method is used for solving, and the initial value of iteration is bestGenerated using a same order explicit method. In actual calculation, in order to avoid increasing calculation time by iteration, iteration is often performed only once, namely, a prediction-correction method is adopted for calculation, and the calculation precision is in the same order as that of the original implicit algorithm. In the invention, a pre-estimation-correction algorithm combining a four-order approximate analysis method and a four-order Gear method is adopted to calculate a re-modeling equation, namely

It has been demonstrated by the former that the order of the Gear method of a-stationary (180 deg. stable region is left half complex plane) does not exceed 2, but the stable region of the four-order Gear method has 146 deg., and the algorithm is also stable as long as the step size is not very large. Simulation results show that the combined algorithm is very effective for solving the rigid ODE initial value problem, matrix inversion is avoided, only one iteration is needed, and the calculation performance is satisfactory.

The combined numerical algorithm solves the initial value problem (19) by the following steps:

□ 1is givenT₀Sufficiently small step size Δ R (which may be taken to be-0.04), and R interval

□ 2 calculating K 'according to the given parameters'_rIs calculated by the formula (23)

And (5) waiting for starting values, and calculating the temperature T by adopting an improved Euler method.

□ 3 recalculated from temperature T'_rAnd is calculated by the formula (26)

Until the R interval is calculated.

Output of

And (5) the final calculation result of T.

(4) Parameter estimation

If the average catalyst activity function phi_iIt is known that the mechanism parameter k in the equation for the reaction rate constant (16) then needs to be estimated before the model equation can be solved_0j，E_j，θ_jFor 17 reactions, this means that there are 51 parameters to be estimated. It is obviously difficult to estimate 51 parameters at the same time, and any algorithm cannot do any job, especially for the case of parameter estimation using industrial data. In order to reduce the parameters to be estimated and the difficulty of parameter estimation, we rewrite the formula (16) into

$k_j=k_{0j}^{\′}\·\exp (-E_j^{\′}/\mathrm{RT})\·P_h^{{\mathrm{\θ}}_j^{\′}}\·\mathrm{\φ},j=1~17----\left(27\right)$

Wherein $k_{0j}^{\′}=k_{0j}\·\exp (-\mathrm{\Δ}{E}_j/\mathrm{RT})\·P_h^{\mathrm{\Δ}{\mathrm{\θ}}_j},$ ΔE_j、Δθ_jRespectively, the estimation error of the corresponding parameter. Activation energy estimation error Δ E for the same catalyst_jGenerally, it does not exceed several tens kJ. mol^-1(ii) a Meanwhile, the reaction pressure drop is generally small (especially a radial reactor), so that k 'is within the reforming reaction temperature range of 450-550℃'_0jSubstantially without large variations with temperature and hydrogen partial pressure. In the present invention, the mechanism model parameters activation energy E and pressure index are both disclosed in the following document [5]]Published data, the parameter to be estimated is only a temperature and pressure independent correction frequency factor k'_0j. Obviously, the number of parameters to be estimated is only 17, thereby greatly reducing the number and difficulty of parameter estimation.

Correcting frequency factor k'_0jThe parameter estimation problem of (2) is actually an unconstrained optimization problem, i.e. the sum of the squares of the errors of the calculated values of the mechanism model and the actual values of the industrial device is required to be minimum. Order to $X={[x_1,x_2\·\·\·,x_n]}^T,x_i=k_{0i}^{\′},$ Where n is 17 as the number of responses, the objective function of the optimization problem can be written as

$\min S\left(X\right)=\frac{1}{P}\underset{k=1}{\overset{p}{\mathrm{\Σ}}}(\underset{i=1}{\overset{q}{\mathrm{\Σ}}}{(y_i^{\mathrm{cal}}-y_i)}^2+\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{(T_i^{\mathrm{cal}}-T_i)}^2)----\left(28\right)$

WhereinFor the calculated values of the model, P is the number of data samples and q is the number of components we are interested in. In the invention, q is 10, namely the total cracked product, hydrogen and P are included₆～P₉₊And A₆～A₉₊And T is the outlet temperature of 4 reactors.

For the solution of equation (28), the Marquardt method (modified Gauss-Newton method) has proven to be simple and effective for the sum of squares form objective functions of the equations. However, this method requires a high initial value and requires the inverse of the matrix, and sometimes fails to converge. The improved variable-scale method (BFGS) not only retains the rapidity of Newton method convergence, but also does not need to ask a second derivative matrix (Hessian matrix) of an objective function and an inverse matrix thereof, and is proved to be one of the most effective unconstrained optimization algorithms. The invention adopts a BFGS method to determine the descending direction of the target function, and then adopts a differential method (one-dimensional search) to obtain the optimal step length.

For the unconstrained optimization problem of equation (28), the BFGS algorithm is as follows:

wherein P is^tThe direction of descent of the objective function s (x) produced for the t-th iteration; f_k(X) is a m ═ q +4 dimensional column vector, the element f of which_k，j(X) (i is 1-P, j is 1-m) is the error of the calculated value and the actual value of each main variable; j. the design is a square_kJacobian matrix being an mxn objective function, J_k＝_F_k(X^(t)) V _X; q is an m x m weight diagonal matrix;

the scale correction matrix used for the t-th iteration is a unit matrix which can be taken as an initial value; lambda [ alpha]_optFor the optimal step size, it can be obtained by an accurate one-dimensional search. Namely have

$S(X^{\left(t\right)}+{\mathrm{\λ}}_{\mathrm{opt}}\·P^t)=\underset{\mathrm{\λ}}{\min }S(X^{\left(t\right)}+\mathrm{\λ}\·P^t)$

Many mature one-dimensional iterative search algorithms (e.g., newton's method, golden section method, polynomial method) can be used to solve this problem. The invention adopts a differential algorithm without iteration, namely S (X)^(t)+λ·P^t) Expanding by Taylor series and taking out the second order term to obtain

$S(X^{\left(t\right)}+\mathrm{\λ}\·P^t)\≈S\left(X^{\left(t\right)}\right)+\mathrm{\ΔS}{\left(X^{\left(t\right)}\right)}^T\·\mathrm{\λ}\·P^t+\frac{1}{2}\mathrm{\λ}\·{\left(P^t\right)}^T\·H(X^{\left(t\right)}\·\mathrm{\λ}\·P^t)----\left(30\right)$

Where H is the Hessian matrix of the objective function. Obviously, lambda is required_optBy making the first derivative of the above formula to λ zero, i.e. by

_S(X^(t))^T·P^t+λ·(P^t)^T·H(X^(t))·P^t＝0 (31)

Then there is

${\mathrm{\λ}}_{\mathrm{opt}}=-\frac{\mathrm{\ΔS}{\left(X^{\left(t\right)}\right)}^T\·P^t}{{\left(P^t\right)}^T\·H\left(X^{\left(t\right)}\right)\·P^t}----\left(32\right)$

Because the objective function in the invention is complex, the Hessian matrix is difficult to solve, but because of the square sum form of the objective function, the product of Jacobian matrix can be used for approximating the Hessian matrix, namely the product of Jacobian matrix

$H\left(X^{\left(t\right)}\right)={\▿}^{2}S\left(X^{\left(t\right)}\right)\≈2\frac{1}{P}\underset{k=1}{\overset{P}{\mathrm{\Σ}}}\{J_k^T\·Q\·J_k\}$

Wherein J_k＝_F_k(X^(t))/_X。

Calculation practice shows that the BFGS adopting the differential accurate one-dimensional search algorithm has higher convergence speed than the BFGS adopting other iterative one-dimensional search algorithms, is not influenced by the initial value of lambda, has the calculation speed equivalent to that of a Marquardt method, avoids solving the inverse of a matrix, greatly reduces the occupation of a memory for high-dimensional problems, and does not have the problem that the matrix is singular and cannot be solved.

Claims (5)

1. A catalytic reforming 17 lumped reaction model modeling method is characterized in that:

(1) subdivision of the reaction feed into 17 lumped components

a) Alkane: methane (C)₁) Ethane (C)₂) Propane (C)₃) N-isomeric butanes (C)₄) N-isomeric Pentanes (C)₅) N-isomeric Hexane (P)₆) N-isomeric Heptane (P)₇) N-iso-octane (P)₈) And an alkane (P) having 9 or more carbon atoms₉₊) The 9 pieces of the total are equal to each other,

b) cycloalkane: methylcyclopentane and cyclohexane (N)₆) Cycloalkane having 7 carbon atoms (N)₇) Cycloalkane having 8 carbon atoms (N)₈) Cycloalkane (N) having 9 or more carbon atoms₉₊) And the 4 pieces of the total are equal,

c) aromatic hydrocarbons: benzene (A)₆) Toluene (A)₇) Xylene and ethylbenzene (A)₈) And heavy aromatic hydrocarbons containing 9 or more carbon atoms (A)₉₊) Equal 4 lumped;

(2) reforming reaction network comprising 17 reactions

The reaction network comprises 4 alkanes (P)₆，P₇，P₈，P₉₊) Dehydrocyclization to give cycloalkane reaction, 4 cycloalkanes (N)₆，N₇，N₈，N₉₊) Dehydrocyclization to generate aromatic hydrocarbon reaction, toluene hydrogenolysis to generate benzene and methane, xylene hydrogenolysis to generate toluene and methane, trimethylbenzene hydrogenolysis to generate xylene and methane, and methylethylbenzene hydrogenolysis to generate toluene and ethane, 4 aromatic hydrocarbon hydrogenolysis reactions, 5 alkane hydrocracking reactions and other 17 reforming reactions;

(3) method for determining model parameters to be estimated

Estimating only 17 corrected frequency factor parameters in the mechanism model;

(4) continuous reforming catalyst coking and deactivation model

The integral term in the catalyst coking model adopts a sectional integral method, and the inactivation model introduces a reactor average activity function, namely the catalyst activity of a single reactor is considered to be the same;

(5) parameter estimation algorithm

A differential accurate one-dimensional search algorithm is introduced on the basis of a classic BFGS algorithm, namely a Hessian matrix of a target function is approximated by utilizing the product of a Jacobian matrix of the target function, so that an analytic expression of the optimal step length of one-dimensional search is solved.

2. The method of claim 1, wherein the reforming reaction network comprising 17 reactions is an alkane dehydrogenation reversible reaction:

P_i_ N_i+H₂，i＝6，7，8，9+naphthene dehydrogenation reversible reaction:

N_i_A_i+3H₂i ═ 6, 7, 8, 9+ alkane hydrocracking reaction:

and (3) carrying out hydrogenolysis reaction of aromatic hydrocarbon:

A₇+H₂→A₆+C₁，A₈+H₂→A₇+C₁

A₉++H₂→A₇+C₂，A₉++H₂→A₈+C₁。

3. the method of claim 1, wherein the 17 modified frequency factor parameters of the model of the estimated-only mechanism are

$k_{0j}^{\′}=k_{0j}\·\exp (-{\mathrm{\ΔE}}_j/\mathrm{RT})\·P_h^{{\mathrm{\Δ\θ}}_j},j=\mathrm{1,2},\·\·\·,17,$

Wherein Δ E_j、Δθ_jRespectively, the estimation error of the corresponding parameter.

4. The method of claim 1, wherein the average activity function of the reactor is

$\stackrel{\‾}{\mathrm{\φ}}\left(t\right)=\underset{N\→\∞}{\lim }\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\φ}}_i/N=\∫\mathrm{\φ}(z,t)\mathrm{dz}/H$

For a radial reactor in any continuous reforming, the catalyst average activity function can be expressed as

${\stackrel{\‾}{\mathrm{\φ}}}_j\left(t\right)={\∫}_{h_j}^{h_j+H_j}\mathrm{\φ}(z,t)\mathrm{dz}/H_j,j=1,\·\·\·,n$

Wherein h is_j＝H_1j-1，H_jRespectively setting the initial value and the whole height of the axial coordinate of the jth reactor, setting the integration time interval as 1h, and the average catalyst activity function expression in each reactor at the time t as

${\stackrel{\‾}{\mathrm{\φ}}}_j\left(t\right)=\frac{1}{t_j}\underset{i=0}{\overset{t_j-1}{\mathrm{\Σ}}}\frac{1}{1+\mathrm{\α}\·\underset{k=0}{\overset{t_{1j-1}+i}{\mathrm{\Σ}}}{\stackrel{\‾}{r}}_c^o(t-k)},t_j=\frac{{\mathrm{\Ω}}_j\·{\mathrm{\ρ}}_c}{v_c}{H}_j,j=1,\·\·\·,n$

It follows that the catalyst average activity function is only a function of the reactor and the run time t, and we consider the catalyst activity function to be the same for 17 reforming reactions, i.e. there is

$r_j=r_j^0\·\stackrel{\‾}{\mathrm{\φ}},0\≤\stackrel{\‾}{\mathrm{\φ}}\≤1,j=\mathrm{1,2},\·\·\·,17.$

5. The method of claim 1, wherein the modeling of the catalytic reforming 17 lumped reaction model is performed in a batch reactor

The analytical formula for determining the optimal step size for the one-dimensional search is

${\mathrm{\λ}}_{\mathrm{opt}}=-\frac{\▿S{\left(X^{\left(t\right)}\right)}^T\·P^t}{{\left(P^t\right)}^T\·H\left(X^{\left(t\right)}\right)\·P^t}$

Wherein the Hessian matrix H (X)^(t)) By products of Jacobian matrices, i.e.

$H\left(X^{\left(t\right)}\right)={\▿}^{2}S\left(X^{\left(t\right)}\right)\≈2\frac{1}{P}\underset{k=1}{\overset{P}{\mathrm{\Σ}}}\{J_k^T\·Q\·J_k\}$

Wherein J_k＝_F_k(X^(t)) X, X are the 17 parameters to be estimated.

Images