CN101776892B - Constraint-prioritized dynamic industrial process optimization system and method - Google Patents

Abstract

The invention relates to a constraint-prioritized dynamic industrial process optimization system, which consists of a field intelligent instrument connected with an industrial process object, a DCS system and an upper computer, wherein the upper computer comprises a constraint conversion module, an initialization module, a double-layered optimization module; the constraint conversion module is used to convert the control variable boundary constraint and the state variable final value constraint in the dynamic fixed boundary optimization problem; the initialization module is used to set initial parameters, discretize a decision vector z(t) and carry out initial assignment; the double-layered optimization module is used to look for decision vectors for optimizing the objective functions of the dynamic optimization problem, and adopts inner-layer optimization and outer-layer optimization for solution, the inner-layer optimization looks for the decision vector z1(t) optimizing the objective function J1, the outer-layer optimization: looking for the decision vector z2(t) optimizing the objective function J2, and an optimal result obtained by the double-layered optimization is stored. The invention also provides a constraint-prioritized dynamic industrial process optimization method. The invention can accurately and quickly find optimal solution for the fixed boundary problem, and has high stability and high applicability.

Description

A kind of constraint-prioritized dynamic industrial process optimization system and method

Technical field

The present invention relates to the optimum control field, especially a kind of constraint-prioritized dynamic industrial process optimization system.

Background technology

Improve constantly along with what the modern industry production run required, dynamic perfromance problems such as non-linear, uncertain, higher-dimension become the bottleneck that enterprise improves output, cuts down the consumption of energy day by day.Adopt dynamic optimization method to solve bottleneck problem and enhancing efficiency by relying on tapping internal latent power in the industrial process optimal control, more and more received the attention of domestic and international academia and industry member.

The dynamic industrial process optimization problem often has state variable boundary values fixed constraint, like the restriction of valve, reactor capacity, pressure, mole fraction etc.Therefore, the boundary values fixation problem is the key and the focus of dynamic industrial process optimization.

Penalty function method is to handle the strategy commonly used of boundary values fixation problem, and it increases the penalty function item and constitutes new objective function on the basis of former objective function, thereby has eliminated the boundary values fixed constraint in the dynamic model; But the value of the validity of penalty function method and penalty factor is closely related; Select not be worthwhilely will cause morbid state, have a strong impact on the calculating effect, do not have established practice and can follow and how to choose suitable penalty factor; Often need progressively tentative calculation, efficient is lower.

Summary of the invention

For overcome traditional penalty function method handle boundary values ill phenomenon can appear during optimization problems fixedly and calculate inaccurate, find the solution inefficient deficiency; The invention provides a kind of optimum solution that can find the boundary values fixation problem accurately and rapidly, and stability is high, applicability is wide constraint-prioritized dynamic industrial process optimization system and method.

The technical solution adopted for the present invention to solve the technical problems is:

A kind of constraint-prioritized dynamic industrial process optimization system; Comprise the field intelligent instrument, DCS system and the host computer that are connected with industrial process object; Described DCS system comprises database and control station; Said field intelligent instrument is connected with said DCS system, and said DCS system is connected with host computer, and described host computer comprises:

The constraint conversion module is used for transforming boundary values fixedly the control variable boundary constraint and the state variable final value constraint (boundary values fixed constraint) of optimization problems, thereby helps finding the solution of optimization problem, takes following steps to accomplish:

(2.1) through intermediate variable processing controls variable edge bound constrained, promptly for having boundary constraint shown in the formula (1):

u _min≤u(t)≤u _max (1)

M ties up control variable u (t), and subscript m in, max represent minimum value and maximal value, u respectively _Min, u _MaxBe constant, the lower bound and the upper bound of the corresponding control variable of difference, take with down conversion:

u(t)＝0.5(u _max-u _min)×{sin[z(t))]+1}+u _min (2)

The trigonometric function that u (t) is converted into the intermediate variable z (t) that does not receive boundary constraint is expressed formula, and finds the solution z (t) as the decision variable of optimization problems;

(2.2) constraint of state variable final value is converted into new objective function, promptly for state variable x with final value constraint _j(t), wherein, c representes to receive the state variable number of final value constraint, x _JfBe given constant, x _j(tf) expression state variable x _j(t) in the value of terminal juncture tf, said final value is constrained to:

x _j(tf)＝x _jf (j＝1，2，...，c) ?(3)

Construct following objective function:

${\mathrm{<figure-callout\; id="1"\; label="J"\; filenames="H2009101566519E00011.png"\; state="\{\{state\}\}">J</figure-callout>}}_1=\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{[x_j\left(\mathrm{tf}\right)-x_{\mathrm{jf}}]}^2---\left(4\right)$

J ₁It is the internal layer objective function that the dual-layer optimization module is found the solution;

Initialization module is used for the setting of initial parameter, discretize and the initial assignment of decision vector z (t), and concrete steps are following:

(3.1) time domain [0, tf] is divided into the N segment: [0, t ₁], [t ₁, t ₂] ..., [t _N-1, t _N], t wherein _N=tf; The length of each time period is tf/N;

(3.2) n dimension decision vector z (t) is carried out discretize on the said time slice of step (3.1), promptly each decision variable is with N the normal value representation of segmentation, and gets initial decision vector z ⁰Be the arbitrary constant vector;

The convergence precision of (3.3) establishing ectonexine optimization is respectively ζ ₁, ζ ₂, when optimization target values iteration error during, stopping iteration less than convergence precision, iterations is respectively k, l, and the initial value of k, l all is taken as 0; If the initial ranging step-length of internal layer optimization is α 0, the initial decision vector of iterative search is z1 ⁰

The dual-layer optimization module, it is optimum to be used to seek the objective function that can not only make optimization problems, with reference to formula (5):

And can satisfy the boundary values fixed constraint, i.e. formula (3) and state equation, i.e. formula (6):

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),z\left(t\right),t],$ x(0)＝x ₀ (6)

Optimizing decision vector z ^*(t), in the formula (5) (6)

ψ is illustrated respectively under the end-condition and at the ingredient of a period of time internal object function, x representes given n dimension state vector, x ₀The expression initial time, i.e. the state vector value of t=0, f representative function vector, the structure of two-layer optimization is found the solution inside and outside taking:

(4.1) internal layer optimization is promptly sought and is made objective function J ₁Optimum decision vector z1 (t), and z1 (t) must satisfy state equation, i.e. the co-state equation of formula (6) and internal layer optimization, i.e. formula (7):

$\frac{\mathrm{d\&lambda;}\left(t\right)}{\mathrm{dt}}=-\mathrm{\&lambda;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;x},$ $\mathrm{\&lambda;}\left(\mathrm{tf}\right)=\frac{{\&PartialD;\mathrm{<figure-callout\; id="1"\; label="J"\; filenames="H2009101566519E00011.png"\; state="\{\{state\}\}">J</figure-callout>}}_1}{\&PartialD;x\left(\mathrm{tf}\right)}---\left(7\right)$

Wherein, λ (t) expression m dimension association attitude vector, subscript T representes vectorial transposition, formula (6) constitutes internal layer ordinary differential equation (Ordinary Differential Equations, ODE) system with formula (7); The optimizing decision variable z1 (t) that internal layer is optimized gained passes to outer as the outer initial solution of optimizing;

(4.2) outer optimization, purpose are on internal layer optimization basis, to search to make objective function J ₂Optimum decision vector z2 (t), and z2 (t) must satisfy state equation, the i.e. co-state equation of formula (6) and outer optimization, i.e. formula (8):

$\frac{\mathrm{d\&theta;}\left(t\right)}{\mathrm{dt}}=\frac{\&PartialD;\mathrm{\&psi;}[x\left(t\right),z\left(t\right),t]}{\&PartialD;x}-\mathrm{\&theta;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;x},$

Wherein, θ (t) expression m dimension association attitude vector,

ψ is illustrated respectively under the end-condition and at a period of time internal object function J ₂Ingredient, formula (6) and formula (8) constitute outer ordinary differential equation ODE system; The outer optimizing decision vector z2 (t) that optimizes gained is exactly the optimum solution z of dual-layer optimization ^*(t), corresponding J ₂Value is exactly the optimal objective value J of dual-layer optimization ^*Then, preserve the optimal result z that dual-layer optimization obtains ^*(t) and J ^*, and it is passed to the output display module.

As preferred a kind of scheme: in the step of said dual-layer optimization module (4.1), internal layer optimization realizes according to following steps:

(4.1.1) choose iteration initial point z1 ⁰, if k=l=0, then z1 ⁰=z ⁰, otherwise z1 ⁰Value is the z2 of outer input ^l

(4.1.2) with the k time iteration point z1 ^k(during k=0, z1 ^k=z1 ⁰) substitution internal layer ODE system, adopt runge kutta method that formula (6) and (7) are carried out forward direction integration and back respectively to integration, solve state vector x and the vectorial λ of association's attitude, and calculate the desired value J of the k time iteration by formula (4) ₁ ^k

(4.1.3) judge the condition of convergence, promptly whether formula (9) is set up, if establishment, then the optimum solution z1 of internal layer ^*=z1 ^k, with z1 ^*Pass to skin, as the initial solution of external iteration; Otherwise change step (4);

|J ₁ ^k-J ₁ ^k+1|≤ζ ₁ (9)

(4.1.4) with state vector x and iteration point z1 ^kSubstitution formula (10) compute gradient g ^k, and preserve z1 ^kAnd g ^k

$g^k\left(t\right)=\mathrm{\&lambda;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}---\left(10\right)$

(4.1.5) confirm best step-size in search α ^k: during as if k=0, get α ^k=α 0, changes step (4.1.6); Otherwise, ask for step factor α according to formula (11)～(13) ^k

${\mathrm{\&alpha;}}^k=\min (\frac{\mathrm{\&pi;}}{D\&CenterDot;\left|\right|g^k\left|\right|},{\mathrm{\&beta;}}^k)---\left(11\right)$

${\mathrm{\&beta;}}^k=\frac{{\left(s^{k-1}\right)}^T\&CenterDot;y^{k-1}}{{\left|\right|y^{k-1}\left|\right|}^2}---\left(12\right)$

s ^k-1＝z1 ^k-z1 ^k-1，y ^k-1＝g ^k-g ^k-1 (13)

Wherein, s ^K-1The error of expression current iteration point and last iteration point, y ^K-1The gradient error of expression current iteration point and last iteration point, D are got the round values in the interval [5,8];

(4.1.6) calculate next iteration point

z1 ^k+1＝z1 ^k-α ^k·g ^k (14)

(4.1.7) iterations is added 1, i.e. k=k+1 is with the z1 in the step (4.1.6) ^K+1Save as current point z1 ^kContinue iteration, change step (4.1.2);

In the said step (4.2), outer optimization according to following algorithm steps realized:

(4.2.1) getting current iteration point is z2 ^l=z1 ^*, l is the current iteration number of times, with the outer ODE of its substitution system, adopts runge kutta method that formula (6) and (8) are carried out forward direction integration and back respectively to integration, solves state vector x and the vectorial θ of association's attitude, and calculates the desired value J of the l time iteration by formula (4) ₂ ^l

(4.2.2) judge the condition of convergence, promptly whether formula (15) is set up, if establishment, then the optimum solution z of dual-layer optimization ^*=z ₂ ^l, optimal objective function value J ^*=J ₂ ^l, preserve and transmit z ^*And J ^*To the output display module; Otherwise change next step;

|J ₂ ^k-J ₂ ^k+1|≤ζ ₂ (15)

(4.2.3) with state vector x and iteration point z1 ^lSubstitution formula (16) compute gradient h ^l, and preserve z2 ^lAnd h ^l

$h^l\left(t\right)=\frac{\&PartialD;\mathrm{\&psi;}[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}+\mathrm{\&theta;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}---\left(16\right)$

(4.2.4) calculate search direction vector d according to formula (17)～(19) ^l, wherein P representes transformation matrix, and subscript l representes iterations, and subscript T representes vectorial transposition, I representation unit matrix, w ^L-1The error of expression current iteration point and last iteration point, v ^L-1The gradient error of expression current iteration point and last iteration point;

d ^l＝-P ^l·h ^l (17)

$P^l=I-\frac{w^{l-1}{\left(v^{l-1}\right)}^T+v^{l-1}{\left(w^{l-1}\right)}^T}{{\left(w^{l-1}\right)}^T{v}^{l-1}}+(1+\frac{{\left|\right|v^{l-1}\left|\right|}^2}{{\left(w^{l-1}\right)}^T{v}^{l-1}})\frac{w^{l-1}{\left(w^{l-1}\right)}^T}{{\left(w^{l-1}\right)}^T{v}^{l-1}}---\left(18\right)$

w ^l-1＝z2 ^l-z2 ^l-1，v ^l-1＝h ^l-h ^l-1 (19)

(4.2.5) carry out precise search, obtain and make J according to Newton method ₂(z2 ^l+ γ d ^l) minimum step-length γ=γ ^l

(4.2.6) calculate next iteration point:

z2 ^l+1＝z2 ^l+γ ^l·d ^l (20)

(4.2.7) iterations is added 1, i.e. l=l+1 is with the z2 in the step (4.2.6) ^L+1Save as current point z2 ^lPassing to internal layer continues to optimize.

Further, said host computer also comprises: information acquisition module, be used to set the sampling time, and gather the multidate information of the industrial process object of uploading by field intelligent instrument.

Further again, said host computer also comprises: the output display module, the optimizing decision that is used for the dual-layer optimization module is calculated is z as a result ^*(t) through type (2) is converted into optimum control path u ^*(t), then with u ^*(t) and optimal objective value J ^*Be transferred to the DCS system, and in the DCS system, show resulting Optimization result information.

A kind of constraint-prioritized dynamic industrial process optimization method, described dynamic optimization method may further comprise the steps:

1) in the DCS system, specifies the state variable and the control variable of dynamic optimization, according to the up-and-down boundary u of the condition enactment control variable of the condition of actual production environment and performance constraint _Max, u _MinWith the sampling period of DCS, and with the historical data of corresponding each variable in the DCS database, control variable up-and-down boundary value u _Max, u _MinSend host computer to;

2) in the constraint conversion module of host computer, take following measure to transform the boundary values fixed constraint in the optimization problems:

(2.1) utilize intermediate variable z (t) to having boundary constraint

u _min≤u(t)≤u _max (1)

M dimension control variable u (t) change:

u(t)＝0.5(u _mam-u _min)×{sin[z(t)]+1}+u _min (2)

The trigonometric function that u (t) is converted into the intermediate variable z (t) that does not receive boundary constraint is expressed formula, and finds the solution z (t) as the decision variable of optimization problems;

(2.2) with the constraint of state variable final value, promptly formula (3) is converted into new objective function J ₁, i.e. formula (4), specific as follows:

x _j(tf)＝x _jf (j＝1，2，...，c) (3)

${\mathrm{<figure-callout\; id="1"\; label="J"\; filenames="H2009101566519E00011.png"\; state="\{\{state\}\}">J</figure-callout>}}_1=\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{[x_j\left(\mathrm{tf}\right)-x_{\mathrm{jf}}]}^2---\left(4\right)$

Wherein, c representes to receive the state variable number of final value constraint, x _JfBe given constant, x _j(tf) expression state variable x _j(t) in the value of terminal juncture tf, J ₁It also is the internal layer objective function that the dual-layer optimization module is found the solution;

3) in the initialization module of host computer, initial parameter is provided with, and the data of DCS system input is carried out initialization process, accomplish according to following steps:

(3.1) time domain [0, tf] is divided into the N segment: [0, t ₁], [t ₁, t ₂] ..., [t _N-1, t _N], t wherein _N=tf; The length of each time period is tf/N;

(3.2) n dimension decision vector z (t) is carried out discretize on (3.1) said time slice, promptly each decision variable is with N the normal value representation of segmentation, and gets initial decision vector z ⁰Be the arbitrary constant vector;

The convergence precision of (3.3) establishing ectonexine optimization is respectively ζ ₁, ζ ₂, iterations is respectively k, l, and the initial value of k, l all is taken as 0; If the initial ranging step-length of internal layer optimization is α 0, the initial decision vector of iterative search is z1 ⁰

4) in the dual-layer optimization module of host computer, searching can not only make the objective function of optimization problems optimum, with reference to formula (5):

And can satisfy the boundary values fixed constraint, i.e. formula (3) and state equation, i.e. formula (6):

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),z\left(t\right),t],$ x(0)＝x ₀ (6)

Optimizing decision vector z ^*(t), in the formula (5) (6)

ψ is illustrated respectively under the end-condition and at the ingredient of a period of time internal object function, x representes given n dimension state vector, x ₀The expression initial time, i.e. the state vector value of t=0, f representative function vector, the structure of two-layer optimization is found the solution inside and outside taking:

(4.1) internal layer optimization is promptly sought and is made objective function J ₁Optimum decision vector z1 (t), and z1 (t) must satisfy state equation, i.e. the co-state equation of formula (6) and internal layer optimization, i.e. formula (7):

$\frac{\mathrm{d\&lambda;}\left(t\right)}{\mathrm{dt}}=-\mathrm{\&lambda;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;x},$ $\mathrm{\&lambda;}\left(\mathrm{tf}\right)=\frac{{\&PartialD;\mathrm{<figure-callout\; id="1"\; label="J"\; filenames="H2009101566519E00011.png"\; state="\{\{state\}\}">J</figure-callout>}}_1}{\&PartialD;x\left(\mathrm{tf}\right)}---\left(7\right)$

Wherein, λ (t) expression m dimension association attitude vector, subscript T representes vectorial transposition, formula (6) constitutes internal layer ordinary differential equation (Ordinary Differential Equations, ODE) system with formula (7); The optimizing decision variable z1 (t) that internal layer is optimized gained passes to outer as the outer initial solution of optimizing;

(4.2) outer optimization, purpose are on internal layer optimization basis, to search to make objective function J ₂Optimum decision vector z2 (t), and z2 (t) must satisfy state equation, the i.e. co-state equation of formula (6) and outer optimization, i.e. formula (8):

$\frac{\mathrm{d\&theta;}\left(t\right)}{\mathrm{dt}}=\frac{\&PartialD;\mathrm{\&psi;}[x\left(t\right),z\left(t\right),t]}{\&PartialD;x}-\mathrm{\&theta;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;x},$

Wherein, θ (t) expression m dimension association attitude vector,

ψ is illustrated respectively under the end-condition and at a period of time internal object function J ₂Ingredient, formula (6) and formula (8) constitute outer ordinary differential equation ODE system; The outer optimizing decision vector z2 (t) that optimizes gained is exactly the optimum solution z of dual-layer optimization ^*(t), corresponding J ₂Value is exactly the optimal objective value J of dual-layer optimization ^*

Then, preserve the optimal result z that dual-layer optimization obtains ^*(t) and J ^*, and it is passed to the output display module.

As preferred a kind of scheme: in the said step (4.1), internal layer optimization realizes according to following steps:

(4.1.1) choose iteration initial point z1 ⁰, if k=l=0, then z1 ⁰=z ⁰, otherwise z1 ⁰Value is the z2 of outer input ^l

(4.1.2) with the k time iteration point z1 ^k(during k=0, z1 ^k=z1 ⁰) substitution internal layer ODE system, adopt runge kutta method that formula (6) and (7) are carried out forward direction integration and back respectively to integration, solve state vector x and the vectorial λ of association's attitude, and calculate the desired value J of the k time iteration by formula (4) ₁ ^k

(4.1.3) judge the condition of convergence, promptly whether formula (9) is set up, if establishment, then the optimum solution z1 of internal layer ^*=z1 ^k, with z1 ^*Pass to skin, as the initial solution of external iteration; Otherwise change step (4);

|J ₁ ^k-J ₁ ^k+1|≤ζ ₁ (9)

(4.1.4) with state vector x and iteration point z1 ^kSubstitution formula (10) compute gradient g ^k, and preserve z1 ^kAnd g ^k

$g^k\left(t\right)=\mathrm{\&lambda;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}---\left(10\right)$

(4.1.5) confirm best step-size in search α ^k: during as if k=0, get α ^k=α 0, changes step (4.1.6); Otherwise, ask for step factor α according to formula (11)～(13) ^k

${\mathrm{\&alpha;}}^k=\min (\frac{\mathrm{\&pi;}}{D\&CenterDot;\left|\right|g^k\left|\right|},{\mathrm{\&beta;}}^k)---\left(11\right)$

${\mathrm{\&beta;}}^k=\frac{{\left(s^{k-1}\right)}^T\&CenterDot;y^{k-1}}{{\left|\right|y^{k-1}\left|\right|}^2}---\left(12\right)$

s ^k-1＝z1 ^k-z1 ^k-1，y ^k-1＝g ^k-g ^k-1 (13)

Wherein, s ^K-1The error of expression current iteration point and last iteration point, y ^K-1The gradient error of expression current iteration point and last iteration point, D are got the round values in the interval [5,8];

(4.1.6) calculate next iteration point

z1 ^k+1＝z1 ^k-α ^k·g ^k (14)

(4.1.7) iterations is added 1, i.e. k=k+1 is with the z1 in the step (4.1.6) ^K+1Save as current point z1 ^kContinue iteration, change step (4.1.2);

In the said step (4.2), outer optimization according to following algorithm steps realized:

(4.2.1) getting current iteration point is z2 ^l=z1 ^*, l is the current iteration number of times, with the outer ODE of its substitution system, adopts runge kutta method that formula (6) and (8) are carried out forward direction integration and back respectively to integration, solves state vector x and the vectorial θ of association's attitude, and calculates the desired value J of the l time iteration by formula (4) ₂ ^l

(4.2.2) judge the condition of convergence, promptly whether formula (15) is set up, if establishment, then the optimum solution z of dual-layer optimization ^*=z2 ^l, optimal objective function value J ^*=J ₂ ^l, preserve and transmit z ^*And J ^*To the output display module; Otherwise change next step;

|J ₂ ^k-J ₂ ^k+1|≤ζ ₂ (15)

(4.2.3) with state vector x and iteration point z1 ^lSubstitution formula (16) compute gradient h ^l, and preserve z2 ^lAnd h ^l

$h^l\left(t\right)=\frac{\&PartialD;\mathrm{\&psi;}[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}+\mathrm{\&theta;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}---\left(16\right)$

(4.2.4) calculate search direction vector d according to formula (17)～(19) ^l, wherein P representes transformation matrix, and subscript l representes iterations, and subscript T representes vectorial transposition, I representation unit matrix, w ^L-1The error of expression current iteration point and last iteration point, v ^L-1The gradient error of expression current iteration point and last iteration point;

d ^l＝-P ^l·h ^l (17)

$P^l=I-\frac{w^{l-1}{\left(v^{l-1}\right)}^T+v^{l-1}{\left(w^{l-1}\right)}^T}{{\left(w^{l-1}\right)}^T{v}^{l-1}}+(1+\frac{{\left|\right|v^{l-1}\left|\right|}^2}{{\left(w^{l-1}\right)}^T{v}^{l-1}})\frac{w^{l-1}{\left(w^{l-1}\right)}^T}{{\left(w^{l-1}\right)}^T{v}^{l-1}}---\left(18\right)$

w ^l-1＝z2 ^l-z2 ^l-1，v ^l-1＝h ^l-h ^l-1 (19)

(4.2.5) carry out precise search, obtain and make J according to Newton method ₂(z2 ^l+ γ d ^l) minimum step-length γ=γ ^l

(4.2.6) calculate next iteration point:

z2 ^l+1＝z2 ^l+γ ^l·d ^l (20)

(4.2.7) iterations is added 1, i.e. l=l+1 is with the z2 in the step (4.2.6) ^L+1Save as current point z2 ^lPassing to internal layer continues to optimize.

Further; Described dynamic optimization method also comprises: the data of the industrial process object that field intelligent instrument is gathered are sent in the real-time data base of DCS system; Output to host computer in each sampling period from the latest data that the database of DCS system obtains, and carry out initialization process at the initialization module of host computer.

Further again, as in said step (4.2.2), to obtain optimizing decision vector z ^*To convert optimum control curve u into through output module as a result ^*(t), and on the man-machine interface of host computer show u ^*(t) and optimal objective value J ^*Simultaneously, optimum control curve u ^*(t) will pass to the control station of DCS system through EBI, and in the DCS system, show resulting Optimization result information.

Technical conceive of the present invention is: the dynamic industrial process optimization problem often has state variable boundary values fixed constraint, and the dynamic optimization method that searching can solve the boundary values fixation problem well is the focus and the difficult point of dynamic industrial process optimization research.Usually the penalty function method that adopts fixedly ill phenomenon possibly occur during optimization problems handling boundary values, has a strong impact on the calculating effect, and all has deficiency aspect accuracy and the counting yield finding the solution.Dynamic optimization method of the present invention; The optimization problems that will have the boundary values fixed constraint through special constraint switching strategy is converted into new dual-layer optimization problem; Combine the algorithm of two kinds of stability and high efficiencies to carry out the iteration optimizing then; When calculating next iteration point, utilized current point with more before any information, have superlinear convergence speed, and can stablize the optimum solution that converges on optimization problems exactly.

Beneficial effect of the present invention mainly shows: the optimum solution that 1, can find the dynamic industrial process optimization problem with boundary values fixed characteristic exactly; 2, efficient height, good stability are found the solution in optimization.Therefore, all be with a wide range of applications in each field of dynamic industrial process optimization.

Description of drawings

Fig. 1 is the hardware structure diagram of dynamic industrial process optimization system provided by the present invention;

Fig. 2 is the functional structure chart that host computer of the present invention is realized dynamic optimization method.

Embodiment

Specify the present invention according to accompanying drawing below.

With reference to Fig. 1, Fig. 2, a kind of constraint-prioritized dynamic industrial process optimization system comprises the field intelligent instrument 2, DCS system and the host computer 6 that are connected with industrial process object 1, and described DCS system is made up of EBI 3, control station 4, database 5; Field intelligent instrument 2 is connected with fieldbus, and said fieldbus is connected with EBI 3, and said EBI is connected with database 5 with control station 4, and described host computer 6 comprises:

Constraint conversion module 8 is used for transforming boundary values fixedly the control variable boundary constraint and the state variable final value constraint (boundary values fixed constraint) of optimization problems, thereby helps finding the solution of optimization problem, takes following steps to accomplish:

(2.1) through intermediate variable processing controls variable edge bound constrained, promptly for having boundary constraint shown in the formula (1)

u _min≤u(t)≤u _max (1)

(subscript m in, max represent minimum value and maximal value, u respectively to m dimension control variable u (t) _Min, u _MaxBe constant, respectively the lower bound and the upper bound of corresponding control variable), take with down conversion:

u(t)＝0.5(u _max-u _min)×{sin[z(t)]+1}+u _min (2)

The trigonometric function that u (t) is converted into the intermediate variable z (t) that does not receive boundary constraint is expressed formula, and finds the solution z (t) as the decision variable of optimization problems;

(2.2) constraint of state variable final value is converted into new objective function, promptly for having the final value constraint

x _j(tf)＝x _jf (j＝1，2，...，c) (3)

State variable x _j(t) (wherein c representes to receive the state variable number of final value constraint, x _JfBe given constant, x _i(tf) expression state variable x _j(t) in the value of terminal juncture tf), construct following objective function:

${\mathrm{<figure-callout\; id="1"\; label="J"\; filenames="H2009101566519E00011.png"\; state="\{\{state\}\}">J</figure-callout>}}_1=\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{[x_j\left(\mathrm{tf}\right)-x_{\mathrm{jf}}]}^2---\left(4\right)$

J ₁The dual-layer optimization module internal layer objective function of finding the solution just.

Initialization module 9 is used for the setting of initial parameter, discretize and the initial assignment of decision vector z (t), and concrete steps are following:

(3.1) time domain [0, tf] is divided into the N segment: [0, t ₁], [t ₁, t ₂] ..., [t _N-1, t _N], t wherein _N=tf; The length of each time period is tf/N;

(3.2) n dimension decision vector z (t) is carried out discretize on (3.1) said time slice, promptly each decision variable is with N the normal value representation of segmentation, and gets initial decision vector z ⁰Be the arbitrary constant vector;

If the convergence precision of ectonexine optimization is respectively ζ ₁, ζ ₂(, stopping iteration) when optimization target values iteration error during less than convergence precision, iterations is respectively k, l (initial value all is taken as 0); If the initial ranging step-length of internal layer optimization is α 0, the initial decision vector of iterative search is z1 ⁰

Dual-layer optimization module 10 is used to seek the objective function that can not only make optimization problems, i.e. formula (5) optimum,

And can satisfy the boundary values fixed constraint, i.e. formula (3) and state equation, i.e. formula (6):

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),z\left(t\right),t],$ x(0)＝x ₀ (6)

Optimizing decision vector z ^*(t), in the formula (5) (6) ψ is illustrated respectively under the end-condition and at the ingredient of a period of time internal object function, x representes given n dimension state vector, x ₀The state vector value of expression initial time (t=0), f representative function vector, the structure of two-layer optimization is found the solution inside and outside taking:

(4.1) internal layer optimization, purpose are to seek to make objective function J ₁Optimum decision vector z1 (t), and z1 (t) must satisfy state equation, and promptly the co-state equation of formula (6) and internal layer optimization is formula (7):

$\frac{\mathrm{d\&lambda;}\left(t\right)}{\mathrm{dt}}=-\mathrm{\&lambda;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;x},$ $\mathrm{\&lambda;}\left(\mathrm{tf}\right)=\frac{{\&PartialD;\mathrm{<figure-callout\; id="1"\; label="J"\; filenames="H2009101566519E00011.png"\; state="\{\{state\}\}">J</figure-callout>}}_1}{\&PartialD;x\left(\mathrm{tf}\right)}---\left(7\right)$

Wherein λ (t) expression m dimension association attitude is vectorial, and subscript T representes vectorial transposition, and formula (6) constitutes internal layer ordinary differential equation (Ordinary Differential Equations, ODE) system with formula (7); The optimizing decision variable z1 (t) that internal layer is optimized gained passes to outer as the outer initial solution of optimizing;

(4.2) outer optimization, purpose are on internal layer optimization basis, to search to make objective function J ₂Optimum decision vector z2 (t), and z2 (t) must satisfy state equation (formula (6)) and the outer co-state equation of optimizing, i.e. formula (8):

$\frac{\mathrm{d\&theta;}\left(t\right)}{\mathrm{dt}}=\frac{\&PartialD;\mathrm{\&psi;}[x\left(t\right),z\left(t\right),t]}{\&PartialD;x}-\mathrm{\&theta;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;x},$

Wherein θ (t) expression m dimension association attitude is vectorial, ψ is illustrated respectively under the end-condition and at a period of time internal object function J ₂Ingredient, formula (6) and formula (8) constitute outer ODE system; The outer optimizing decision vector z2 (t) that optimizes gained is exactly the optimum solution z of dual-layer optimization ^*(t), corresponding J ₂Value is exactly the optimal objective value J of dual-layer optimization ^*

Then, preserve the optimal result z that dual-layer optimization obtains ^*(t) and J ^*, and with its pass to output display module 11.

The dual-layer optimization module 10 of said host computer 6 adopts different algorithms to carry out ectonexine optimization respectively.Internal layer optimization realizes according to following algorithm steps:

(4.1.1) choose iteration initial point z1 ⁰, if k=l=0, then z1 ⁰=z ⁰, otherwise z1 ⁰Value is the z2 of outer input ^l

(4.1.2) with the k time iteration point z1 ^k(during k=0, z1 ^k=z1 ⁰) substitution internal layer ODE system, adopt runge kutta method that formula (6) and (7) are carried out forward direction integration and back respectively to integration, solve state vector x and the vectorial λ of association's attitude, and calculate the desired value J of the k time iteration by formula (4) ₁ ^k

(4.1.3) judge the condition of convergence, promptly whether formula (9) is set up, if establishment, then the optimum solution z1 of internal layer ^*=z1 ^k, with z1 ^*Pass to skin, as the initial solution of external iteration; Otherwise change step (4);

|J ₁ ^k-J ₁ ^k+1|≤ζ ₁ (9)

(4.1.4) with state vector x and iteration point z1 ^kSubstitution formula (10) compute gradient g ^k, and preserve z1 ^kAnd g ^k

$g^k\left(t\right)=\mathrm{\&lambda;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;z}---\left(10\right)$

(4.1.5) confirm best step-size in search α ^k: during as if k=0, get α ^k=α 0, changes step (6); Otherwise, ask for step factor α according to formula (11)～(13) ^k

${\mathrm{\&alpha;}}^k=\min (\frac{\mathrm{\&pi;}}{D\&CenterDot;\left|\right|g^k\left|\right|},{\mathrm{\&beta;}}^k)---\left(11\right)$

${\mathrm{\&beta;}}^k=\frac{{\left(s^{k-1}\right)}^T\&CenterDot;y^{k-1}}{{\left|\right|y^{k-1}\left|\right|}^2}---\left(12\right)$

s ^k-1＝z1 ^k-z1 ^k-1，y ^k-1＝g ^k-g ^k-1 (13)

Wherein, s ^K-1The error of expression current iteration point and last iteration point, y ^K-1Represent current

The gradient error of iteration point and last iteration point, D are got the round values in the interval [5,8];

(4.1.6) calculate next iteration point

z1 ^k+1＝z1 ^k-α ^k·g ^k (14)

(4.1.7) iterations is added 1, i.e. k=k+1 is with the z1 in the step (4.1.6) ^K+1Save as current point z1 ^kContinue iteration, change step (4.1.2).

Outer optimization according to following algorithm steps realized:

(4.2.1) getting current iteration point is z2 ^l=z1 ^*(l is the current iteration number of times) with the outer ODE of its substitution system, adopts runge kutta method that formula (6) and (8) are carried out forward direction integration and back respectively to integration, solves state vector x and the vectorial θ of association's attitude, and calculates the desired value J of the l time iteration by formula (4) ₂ ^l

(4.2.2) judge the condition of convergence, promptly whether formula (15) is set up, if establishment, then the optimum solution z of dual-layer optimization ^*=z2 ^l, optimal objective function value J ^*=J ₂ ^l, preserve and transmit z ^*And J ^*To the output display module; Otherwise change next step;

|J ₂ ^k-J ₂ ^k+1|≤ζ ₂ (15)

(4.2.3) with state vector x and iteration point z1 ^lSubstitution formula (16) compute gradient h ^l, and preserve z2 ^lAnd h ^l

$h^l\left(t\right)=\frac{\&PartialD;\mathrm{\&psi;}[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}+\mathrm{\&theta;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}---\left(16\right)$

(4.2.4) calculate search direction vector d according to formula (17)～(19) ^l, wherein P representes transformation matrix, and subscript l representes iterations, and subscript T representes vectorial transposition, I representation unit matrix, w ^L-1The error of expression current iteration point and last iteration point, v ^L-1The gradient error of expression current iteration point and last iteration point;

d ^l＝-P ^l·h ^l (17)

$P^l=I-\frac{w^{l-1}{\left(v^{l-1}\right)}^T+v^{l-1}{\left(w^{l-1}\right)}^T}{{\left(w^{l-1}\right)}^T{v}^{l-1}}+(1+\frac{{\left|\right|v^{l-1}\left|\right|}^2}{{\left(w^{l-1}\right)}^T{v}^{l-1}})\frac{w^{l-1}{\left(w^{l-1}\right)}^T}{{\left(w^{l-1}\right)}^T{v}^{l-1}}---\left(18\right)$

w ^l-1＝z2 ^l-z2 ^l-1，v ^l-1＝h ^l-h ^l-1 (19)

(4.2.5) carry out precise search, obtain and make J according to Newton method ₂(z2 ^l+ γ d ^l) minimum step-length γ=γ ^l

(4.2.6) calculate next iteration point:

z2 ^l+1＝z2 ^l+γ ^l·d ^l (20)

(4.2.7) iterations is added 1, i.e. l=l+1 is with the z2 in the step (4.2.6) ^L+1Save as current point z2 ^lPassing to internal layer continues to optimize.

Said host computer 6 also comprises: information acquisition module 7, be used to set the sampling time, and gather the multidate information of the industrial process object of uploading by field intelligent instrument.

Said host computer 6 also comprises: output display module 11, the optimizing decision that is used for dual-layer optimization module 10 is calculated is z as a result ^*(t) through type (2) is converted into optimum control path u ^*(t), then with u ^*(t) and optimal objective value J ^*Be transferred to the DCS system, and in the DCS system, show resulting Optimization result information.

The system hardware structure figure of this case study on implementation is shown in accompanying drawing 1; The core of said dynamic optimization system comprises constraint conversion module 8, initialization module 9 and the dual-layer optimization module 10 3 big functional modules in the host computer 6 of being with man-machine interface, comprises in addition: field intelligent instrument 2, DCS system and fieldbus.Described DCS system is made up of EBI 3, control station 4, database 5; Industrial process object 1, field intelligent instrument 2, DCS system, host computer 6 link to each other through fieldbus successively, realize uploading and assigning of information flow, and host computer and first floor system are in time carried out message exchange, realize the on-line optimization of system.

Embodiment 2

See figures.1.and.2, a kind of constraint-prioritized dynamic industrial process optimization method, described dynamic optimization method is implemented according to following steps:

1), in the DCS system, specifies the state variable and the control variable of dynamic optimization, according to the up-and-down boundary u of the condition enactment control variable of the condition of actual production environment and performance constraint _Max, u _MinWith the sampling period of DCS, and with the historical data of corresponding each variable in the DCS database 5, control variable up-and-down boundary value u _Max, u _MinSend host computer 6 to.

2), in the constraint conversion module 8 of host computer, at first, utilize intermediate variable z (t) that the m dimension control variable u (t) with boundary constraint shown in the formula (1) is changed, i.e. formula (2),

u _min≤u(t)≤u _max (1)

u(t)＝0.5(u _max-u _min)×{sin[z(t)]+1}+u _min (2)

The trigonometric function that u (t) is converted into the intermediate variable z (t) that does not receive boundary constraint is expressed formula, and finds the solution z (t) as the decision variable of optimization problems; Then, with the constraint of state variable final value, promptly formula (3) is converted into new objective function J ₁, i.e. formula (4):

x _j(tf)＝x _if (j＝1，2，...，c) (3)

${\mathrm{<figure-callout\; id="1"\; label="J"\; filenames="H2009101566519E00011.png"\; state="\{\{state\}\}">J</figure-callout>}}_1=\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{[x_j\left(\mathrm{tf}\right)-x_{\mathrm{jf}}]}^2---\left(4\right)$

Wherein, c representes to receive the state variable number of final value constraint, x _JfBe given constant, x _j(tf) expression state variable x _j(t) in the value of terminal juncture tf, J ₁It also is the internal layer objective function that dual-layer optimization module 10 is found the solution.

3) in the initialization module 9 of host computer, initial parameter is provided with, and the data of DCS system input is carried out initialization process, accomplish according to following steps:

(3.1) time domain [0, tf] is divided into the N segment: [0, t ₁], [t ₁, t ₂] ..., [t _N-1, t _N], t wherein _N=tf; The length of each time period is tf/N;

(3.2) n dimension decision vector z (t) is carried out discretize on (3.1) said time slice, promptly each decision variable is with N the normal value representation of segmentation, and gets initial decision vector z ⁰Be the arbitrary constant vector;

The convergence precision of (3.3) establishing ectonexine optimization is respectively ζ ₁, ζ ₂, iterations is respectively k, l, and k, l initial value all are taken as 0; If the initial ranging step-length of internal layer optimization is α 0, the initial decision vector of iterative search is z1 ⁰

4) in the dual-layer optimization module 10 of host computer, searching can not only make the objective function J of optimization problems ₂, i.e. formula (5) optimum, and can satisfy the boundary values fixed constraint, i.e. formula (3) and state equation, the i.e. optimizing decision of formula (6) vector z ^*(t), and with z ^*(t) and corresponding optimal objective value J ^*Pass to the output display module; Structure through two-layer optimization inside and outside taking is found the solution;

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),z\left(t\right),t],$ x(0)＝x ₀ (6)

(4.1) internal layer optimization, purpose are to seek to make objective function J ₁Optimum decision vector z1 (t), and z1 (t) must satisfy the co-state equation (formula (7)) of state equation (formula (6)) and internal layer optimization:

$\frac{\mathrm{d\&lambda;}\left(t\right)}{\mathrm{dt}}=-\mathrm{\&lambda;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;x},$ $\mathrm{\&lambda;}\left(\mathrm{tf}\right)=\frac{{\&PartialD;\mathrm{<figure-callout\; id="1"\; label="J"\; filenames="H2009101566519E00011.png"\; state="\{\{state\}\}">J</figure-callout>}}_1}{\&PartialD;x\left(\mathrm{tf}\right)}---\left(7\right)$

(4.2) outer optimization, purpose are on internal layer optimization basis, to search to make objective function J ₂Optimum decision vector z2 (t), and z2 (t) must satisfy state equation (formula (6)) and the outer co-state equation of optimizing (formula (8)):

$\frac{\mathrm{d\&theta;}\left(t\right)}{\mathrm{dt}}=\frac{\&PartialD;\mathrm{\&psi;}[x\left(t\right),z\left(t\right),t]}{\&PartialD;x}-\mathrm{\&theta;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;x},$

In the said step (4.1), internal layer optimization realizes according to following steps:

(4.1.1) choose iteration initial point z1 ⁰, if k=l=0, then z1 ⁰=z ⁰, otherwise z1 ⁰Value is the z2 of outer input ^l

(4.1.2) with the k time iteration point z1 ^k(during k=0, z1 ^k=z1 ⁰) substitution internal layer ODE system, adopt runge kutta method that formula (6) and (7) are carried out forward direction integration and back respectively to integration, solve state vector x and the vectorial λ of association's attitude, and calculate the desired value J of the k time iteration by formula (4) ₁ ^k

(4.1.3) judge the condition of convergence, promptly whether formula (9) is set up, if establishment, then the optimum solution z1 of internal layer ^*=z1 ^k, with z1 ^*Pass to skin, as the initial solution of external iteration; Otherwise change step (4);

|J ₁ ^k-J ₁ ^k+1|≤ζ ₁ (9)

(4.1.4) with state vector x and iteration point z1 ^kSubstitution formula (10) compute gradient g ^k, and preserve z1 ^kAnd g ^k

$g^k\left(t\right)=\mathrm{\&lambda;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}---\left(10\right)$

(4.1.5) confirm best step-size in search α ^k: during as if k=0, get α ^k=α 0, changes step (4.1.6); Otherwise, ask for step factor α according to formula (11)～(13) ^k

${\mathrm{\&alpha;}}^k=\min (\frac{\mathrm{\&pi;}}{D\&CenterDot;\left|\right|g^k\left|\right|},{\mathrm{\&beta;}}^k)---\left(11\right)$

${\mathrm{\&beta;}}^k=\frac{{\left(s^{k-1}\right)}^T\&CenterDot;y^{k-1}}{{\left|\right|y^{k-1}\left|\right|}^2}---\left(12\right)$

s ^k-1＝z1 ^k-z1 ^k-1，y ^k-1＝g ^k-g ^k-1 (13)

Wherein, s ^K-1The error of expression current iteration point and last iteration point, y ^K-1Represent current

The gradient error of iteration point and last iteration point, D are got the round values in the interval [5,8];

(4.1.6) calculate next iteration point

z1 ^k+1＝z1 ^k-α ^k·g ^k (14)

(4.1.7) iterations is added 1, i.e. k=k+1 is with the z1 in the step (4.1.6) ^K+1Save as current point z1 ^kContinue iteration, change step (4.1.2);

In the said step (4.2), outer optimization according to following algorithm steps realized:

(4.2.1) getting current iteration point is z2 ^l=z1 ^*, l is the current iteration number of times, with the outer ODE of its substitution system, adopts runge kutta method that formula (6) and (8) are carried out forward direction integration and back respectively to integration, solves state vector x and the vectorial θ of association's attitude, and calculates the desired value J of the l time iteration by formula (4) ₂ ^l

(4.2.2) judge the condition of convergence, promptly whether formula (15) is set up, if establishment, then the optimum solution z of dual-layer optimization ^*=z2 ^l, optimal objective function value J ^*=J ₂ ^l, preserve and transmit z ^*And J ^*To the output display module; Otherwise change next step;

|J ₂ ^k-J ₂ ^k+1|≤ζ ₂ (15)

(4.2.3) with state vector x and iteration point z1 ^lSubstitution formula (16) compute gradient h ^l, and preserve z2 ^lAnd h ^l

$h^l\left(t\right)=\frac{\&PartialD;\mathrm{\&psi;}[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}+\mathrm{\&theta;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}---\left(16\right)$

(4.2.4) calculate search direction vector d according to formula (17)～(19) ^l, wherein P representes transformation matrix, and subscript l representes iterations, and subscript T representes vectorial transposition, I representation unit matrix, w ^L-1The error of expression current iteration point and last iteration point, v ^L-1The gradient error of expression current iteration point and last iteration point;

d ^l＝-P ^l.h ^l (17)

$P^l=I-\frac{w^{l-1}{\left(v^{l-1}\right)}^T+v^{l-1}{\left(w^{l-1}\right)}^T}{{\left(w^{l-1}\right)}^T{v}^{l-1}}+(1+\frac{{\left|\right|v^{l-1}\left|\right|}^2}{{\left(w^{l-1}\right)}^T{v}^{l-1}})\frac{w^{l-1}{\left(w^{l-1}\right)}^T}{{\left(w^{l-1}\right)}^T{v}^{l-1}}---\left(18\right)$

w ^l-1＝z2 ^l-z2 ^l-1，v ^l-1＝h ^l-h ^l-1 (19)

(4.2.5) carry out precise search, obtain and make J according to Newton method ₂(z2 ^l+ γ d ^l) minimum step-length

γ＝γ ^l；

(4.2.6) calculate next iteration point:

z2 ^l+1＝z2 ^l+γ ^l·d ^l (20)

(4.2.7) iterations is added 1, i.e. l=l+1 is with the z2 in the step (4.2.6) ^L+1Save as current point z2 ^lPassing to internal layer continues to optimize.

In the information acquisition module 7 of host computer; Set the sampling time; The data of the industrial process object of collection site intelligence instrument input; Be saved in the real-time data base 5 of DCS system, and the constraint conversion module 8 that the latest data of database 5 gained of DCS system outputs to host computer handled with initialization module 9 in each sampling period;

In the output display module 11 of host computer, the optimizing decision of gained vector z in the said step (4-B2) ^*Convert through type (2) into optimum control curve u ^*(t), and with optimal objective value J ^*Be presented at together on the man-machine interface of host computer; Simultaneously, optimum control curve u ^*(t) will pass to the control station 4 of DCS system through EBI, and in the DCS system, show resulting Optimization result information.

In the present embodiment, system begins to put into operation, and detailed process is:

1) utilizes timer, set the time interval of each Data Detection and collection;

2) field intelligent instrument 2 detects the data of industrial process object 1 and is sent to the DCS system

In the real-time data base 5, obtain up-to-date variable data;

3) in the constraint conversion module 8 of host computer 6, boundary constraint is handled to control variable,

With the input of process result as initialization module 9 and dual-layer optimization module 10;

4) in the initialization module 9 of host computer 6, according to actual production demand and performance constraint condition each module correlation parameter and variable are carried out initialization process, with the input of process result as dual-layer optimization module 10;

5) the dual-layer optimization module 10 of host computer 6 is carried out double-deck iteration optimization according to the substitution of variable relation and the new objective function information of constraint conversion module 8, and results of optimization is sent to output display module 11;

6) the output display module 11 of host computer 6; Substitution of variable relation according to constraint conversion module 8 is changed the optimizing decision curve that dual-layer optimization module 10 calculates; Optimum control object information with gained is transferred to the DCS system then; And be shown in the control station 4 of the man-machine interface and the DCS system of host computer 6, through DCS system and fieldbus resulting Optimization result information transmission is shown to the work on the spot station simultaneously, and carry out optimum operation by the work on the spot station.

The foregoing description is used for the present invention that explains, rather than limits the invention, and in the protection domain of spirit of the present invention and claim, any modification and change to the present invention makes all fall into protection scope of the present invention.

Claims (3)

1. constraint-prioritized dynamic industrial process optimization system; Comprise the field intelligent instrument, DCS system and the host computer that are connected with industrial process object; Described DCS system comprises database and control station; Said field intelligent instrument is connected with said DCS system, and said DCS system is connected with host computer, and it is characterized in that: described host computer comprises:

The constraint conversion module is used for transforming boundary values fixedly the control variable boundary constraint and the constraint of state variable final value of optimization problems, takes following steps to accomplish:

(2.1) through intermediate variable processing controls variable edge bound constrained, promptly for having boundary constraint shown in the formula (1):

u _min≤u(t)≤u _max (1)

M ties up control variable u (t), and subscript m in, max represent minimum value and maximal value, u respectively _Min, u _MaxBe constant, the lower bound and the upper bound of the corresponding control variable of difference, take with down conversion:

u(t)＝0.5(u _max-u _min)×{sin[z(t)]+1}+u _min (2)

The trigonometric function that u (t) is converted into the intermediate variable z (t) that does not receive boundary constraint is expressed formula, and finds the solution z (t) as the decision vector of optimization problems;

(2.2) constraint of state variable final value is converted into new objective function, promptly for state variable x with final value constraint _j(t), construct following objective function:

$J_1=\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{[x_j\left(\mathrm{tf}\right)-x_{\mathrm{jf}}]}^2---\left(3\right)$

Wherein, J ₁Be the internal layer objective function that the dual-layer optimization module is found the solution, c representes to receive the state variable number of final value constraint, x _JfBe given constant, x _j(tf) expression state variable x _j(t) in the value of terminal juncture tf, said final value is constrained to:

x _j(tf)＝x _jf(j＝1，2，....，c) (4)

Initialization module is used for the setting of initial parameter, discretize and the initial assignment of decision vector z (t), and concrete steps are following:

(3.1) time domain [0, tf] is divided into the N segment: [0, t ₁], [t ₁, t ₂] ..., [t _N-1, t _N], t wherein _N=tf; The length of each time period is tf/N;

(3.2) n dimension decision vector z (t) is carried out discretize on the said time slice of step (3.1), promptly each decision vector is with N the normal value representation of segmentation, and gets initial decision vector z ⁰Be the arbitrary constant vector;

The convergence precision of (3.3) establishing ectonexine optimization is respectively ζ ₁, ζ ₂, when optimization target values iteration error during, stopping iteration less than convergence precision, iterations is respectively k, q, and the initial value of k, q all is taken as 0; If the initial ranging step-length of internal layer optimization is α 0, the initial decision vector of iterative search is z1 ⁰

The dual-layer optimization module, it is optimum to be used to seek the objective function that can not only make optimization problems, with reference to formula (5):

And can satisfy final value constraint, i.e. formula (3) and state equation, i.e. formula (6):

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),z\left(t\right),t],x\left(0\right)=x_0---\left(6\right)$

Optimizing decision vector z ^*(t), in the formula (5) (6)

ψ is illustrated respectively under the end-condition and at the ingredient of a period of time internal object function, x representes given n dimension state vector, x ₀The expression initial time, i.e. the state vector value of t=0, f representative function vector, the structure of two-layer optimization is found the solution inside and outside taking:

(4.1) internal layer optimization is promptly sought and is made objective function J ₁Optimum internal layer optimizing decision vector z1 (t), and z1 (t) must satisfy state equation, i.e. the co-state equation of formula (6) and internal layer optimization, i.e. formula (7):

$\frac{\mathrm{d\&lambda;}\left(t\right)}{\mathrm{dt}}=-\mathrm{\&lambda;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;x},\mathrm{\&lambda;}\left(\mathrm{tf}\right)=\frac{\&PartialD;J_1}{\&PartialD;x\left(\mathrm{tf}\right)}---\left(7\right)$

Wherein, λ (t) expression m dimension association attitude vector, subscript T representes vectorial transposition, formula (6) constitutes the internal layer coefficient ordinary differential equation system with formula (7); The internal layer optimizing decision vector z1 (t) that internal layer is optimized gained passes to outer as the outer initial solution of optimizing;

(4.2) outer optimization, purpose are on internal layer optimization basis, to search to make objective function J ₂Optimum outer optimizing decision vector z2 (t), and z2 (t) must satisfy state equation, the i.e. co-state equation of formula (6) and outer optimization, i.e. formula (8):

Wherein, θ (t) expression m dimension association attitude vector,

ψ is illustrated respectively under the end-condition and at a period of time internal object function J ₂Ingredient, formula (6) and formula (8) constitute outer coefficient ordinary differential equation system; The optimizing decision vector z that the outer outer optimizing decision vector z2 (t) that optimizes gained is exactly a dual-layer optimization ^*(t), corresponding J ₂Value is exactly the optimal objective value J of dual-layer optimization ^*

Then, preserve the optimal result z that dual-layer optimization obtains ^*(t) and J ^*, and it is passed to the output display module;

The output display module is used for the optimizing decision vector z that the dual-layer optimization module is calculated ^*(t) through type (2) is converted into optimum control path u ^*(t), then with u ^*(t) and optimal objective value J ^*Be transferred to the DCS system, and in the DCS system, show resulting Optimization result information.

2. a kind of constraint-prioritized dynamic industrial process optimization system according to claim 1 is characterized in that: in the said step (4.1), internal layer optimization realizes according to following steps:

(4.1.1) choose iteration initial point z1 ⁰, if k=q=0, then z1 ⁰=z ⁰, otherwise z1 ⁰Value is the z2 of outer input ^q

(4.1.2) with the k time iteration point z1 ^kSubstitution internal layer ODE system, during k=0, z1 ^k=z1 ⁰, adopt runge kutta method that formula (6) and (7) are carried out forward direction integration and back respectively to integration, solve state vector x and the vectorial λ of association's attitude, and calculate the desired value J of the k time iteration by formula (4) ₁ ^k

(4.1.3) judge the condition of convergence, promptly whether formula (9) is set up, if set up, then internal layer optimizing decision vector separates z1 ^*(t)=z1 ^k, with z1 ^*(t) pass to skin, as the initial solution of external iteration; Otherwise change step (4.1.4);

|J ₁ ^k-J ₁ ^k+1|＜ζ ₁ (9)

(4.1.4) with state vector x and iteration point z1 ^kSubstitution formula (10) compute gradient g ^k, and preserve z1 ^kAnd g ^k

$g^k\left(t\right)=\mathrm{\&lambda;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}---\left(10\right)$

(4.1.5) confirm best step-size in search α ^k: during as if k=0, get α ^k=α 0, changes step (4.1.6); Otherwise, ask for step factor α according to formula (11)～(13) ^k

${\mathrm{\&alpha;}}^k=\min (\frac{\mathrm{\&pi;}}{D\&CenterDot;\left|\right|g^k\left|\right|},{\mathrm{\&beta;}}^k)---\left(11\right)$

${\mathrm{\&beta;}}^k=\frac{{\left(s^{k-1}\right)}^T\&CenterDot;y^{k-1}}{{\left|\right|y^{k-1}\left|\right|}^2}---\left(12\right)$

s ^k-1＝z1 ^k-z1 ^k-1，y ^k-1＝g ^k-g ^k-1 (13)

Wherein, s ^K-1The error of expression current iteration point and last iteration point, y ^K-1The gradient error of expression current iteration point and last iteration point, D are got the round values in the interval [5,8];

(4.1.6) calculate next iteration point

z1 ^k+1＝z1 ^k-α ^k·g ^k (14)

(4.1.7) iterations is added 1, i.e. k=k+1 is with the z1 in the step (4.1.6) ^K+1Save as current point z1 ^kContinue iteration, change step (4.1.2);

In the said step (4.2), outer optimization according to following algorithm steps realized:

(4.2.1) getting current iteration point is z2 ^q=z1 ^*, q is the current iteration number of times, with the outer ODE of its substitution system, adopts runge kutta method that formula (6) and (8) are carried out forward direction integration and back respectively to integration, solves state vector x and the vectorial θ of association's attitude, and calculates the desired value J of the q time iteration by formula (4) ₂ ^q

(4.2.2) judge the condition of convergence, promptly whether formula (15) is set up, if set up, then the optimizing decision of dual-layer optimization vector separates z ^*(t)=z2 ^q, optimal objective function value J ^*=J ₂ ^q, preserve and transmit z ^*(t) and J ^*To the output display module; Otherwise change next step;

|J ₂ ^q-J ₂ ^q+1|＜ζ ₂ (15)

(4.2.3) with state vector x and iteration point z2 ^qSubstitution formula (16) compute gradient h ^q, and preserve z2 ^qAnd h ^q

$h^q\left(t\right)=\frac{\&PartialD;\mathrm{\&psi;}[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}+\mathrm{\&theta;}{\left(t\right)}^T\&CenterDot;\frac{\&PartialD;f[x\left(t\right),z\left(t\right),t]}{\&PartialD;z\left(t\right)}---\left(16\right)$

(4.2.4) calculate search direction vector d according to formula (17)～(19) ^q, wherein P representes transformation matrix, and subscript q representes iterations, and subscript T representes vectorial transposition, I representation unit matrix, w ^Q-1The error of expression current iteration point and last iteration point, v ^Q-1The gradient error of expression current iteration point and last iteration point;

d ^q＝-P ^q.h ^q (17)

$P^q=I-\frac{w^{q-1}{\left(v^{q-1}\right)}^T+v^{q-1}{\left(w^{q-1}\right)}^T}{{\left(w^{q-1}\right)}^T{v}^{q-1}}+(1+\frac{{\left|\right|v^{q-1}\left|\right|}^2}{{\left(w^{q-1}\right)}^T{v}^{q-1}})\frac{w^{q-1}{\left(w^{q-1}\right)}^T}{{\left(w^{q-1}\right)}^T{v}^{q-1}}---\left(18\right)$

w ^q-1＝z2 ^q-z2 ^q-1，v ^q-1＝h ^q-h ^q-1

(19)

(4.2.5) carry out precise search, obtain and make J according to Newton method ₂(z2 ^q+ γ d ^q) minimum step-length γ=γ ^q

(4.2.6) calculate next iteration point:

z2 ^q+1＝z2 ^q+γ ^q·d ^q (20)

(4.2.7) iterations is added 1, i.e. q=q+1 is with the z2 in the step (4.2.6) ^Q+1Save as current point z2 ^qPassing to internal layer continues to optimize.

3. according to claim 1 or claim 2 a kind of constraint-prioritized dynamic industrial process optimization system; It is characterized in that: said host computer also comprises: information acquisition module; Be used to set the sampling time, gather the multidate information of the industrial process object of uploading by field intelligent instrument.

Images