CN104281076B - Control method of protein concentration - Google Patents

Abstract

A controlling method for the concentration of protein includes the steps that (1), system identification is carried out on a controlled object through related chemical reactions in a system to obtain a transfer function model of the controlled object; (2), values of time delay coefficients tau1, tau2 and tau3 of a controller are selected appropriately; (3), the characteristic function of the closed-loop system is solved; (4), the value of k1 is fixed; (5), an appropriate weight function W(s) is selected, the compensation sensitivity function T (s, k1, k2 and k3) of the closed-loop system is calculated, and a suitable H8 performance index gamma is selected, so that W(s) and T (s, k1, k2 and k3) meet an expression (please see the expression in the specification), the characteristic function corresponding to the expression is solved, the boundary line, on the plane (k2 and k3), of the H8 performance index is solved, which side the parameter domain of the H8 performance index is located on is judged, the intersection set between the region judged to meet the H8 performance index and a stability region is gotten, and therefore the intersection set region is a control parameter region meeting the H8 performance index.

Description

Control method of protein concentration

technical field

The invention relates to a method for controlling protein concentration, in particular to an _H∞ distributed time-delay control method for a protein synthesis system with time-delay.

Background technique

Protein is the scaffold and main substance that constitutes the tissues and organs of the human body, and plays a very important role in the life activities of the human body. It can be said that there is no life activity without protein. As the bearer of life activities, the realization of most of the physiological functions of organisms depends on protein. Protein is the carrier of biological functions, and has the functions of catalysis (such as enzyme), regulation (such as insulin), immunity (such as antibody), movement ( Such as myosin and actin in muscle fibers), transport (such as hemoglobin), control (such as repressor proteins), structural components (such as structural proteins), storage, scaffolding, defense and offense and other functions. Both protein deficiency and excessive protein intake will cause damage to the living body, and even cause certain bodily functions of the living body to be disordered.

Currently, protein synthesis systems are mainly used in recombinant biological cells or to produce polypeptides or proteins of interest in cells. The present invention relates to the design of _H∞ distributed time-delay controller of protein synthesis system with time-delay. Murry, an expert in the field of control, provides a method for the design and simulation of protein concentration regulation systems, and controls the synthesis of target proteins by using a single negative feedback loop, so that the protein concentration reaches a set concentration (de los Santos ELC, HsiaoV, Murray R M. DESIGN AND IMPLEMENTATION OF A BIOMOLECULAR CIRCUIT FOR TRACKING PROTEIN CONCENTRATION[J]). Although this method constitutes a closed-loop system, it does not add a controller, so that the system performance has not achieved the effects of small overshoot, rapid stabilization, and strong anti-interference. Therefore, it is necessary to introduce a distributed time-delay controller to control the protein concentration. Distributed time-delay controllers actively introduce time-delays, which are generally considered to be unfavorable factors, to improve the stability and performance of control systems. The advantage of introducing time-delay in the controller is to convert the time-delay existing in the system into favorable factors. It has the advantages of simple mechanism, low computational complexity, and easy implementation. It can be used in occasions where ordinary controllers are difficult to implement (such as gene regulation networks). , realize the asynchronous operation of the system, easily transform into a digital controller, realize the direct control of the continuous system, and have the characteristics similar to the hybrid system.

A stabilization system is a prerequisite for designing a control system. Singh proposed the use of single time-delay and distributed time-delay controllers to control the system, and gave a method for tuning parameters (Singh T, Vadali S R. Robust time-delay control [J]. ASME Journal of Dynamic Systems, Measurement , and Control, 1993, 115(2A):303-306). However, its parameter setting method is only designed for the second-order oscillation link, which has great limitations, because in actual industrial control, the actual system does not necessarily only contain the second-order oscillation link, and has time-delay characteristics. Gu provides a solution method with geometric characteristics and detailed parameters of multi-delay linear systems in the literature (Gu K, Naghnaeian M. Stability crossing set for systems with three delays[J]. Automatic Control, IEEE Transactions on, 2011, 56(1):11-26). Although the parameter domain that makes the system stable is solved, the performance of the system is not required. How to design a distributed time-delay controller in a feedback control system to ensure the stability of the time-delay controlled system and make the system performance meet the H _∞ performance index is an unsolved research problem that plays an important role in the field of control.

Contents of the invention

The present invention overcomes the above-mentioned deficiencies of the prior art, and provides a protein concentration control method that stabilizes the distributed time-delay control of a protein production system with a time-delay.

The present invention first uses the improved D-segmentation method and the determination method of which side of the corresponding boundary line is the stable region to determine the stable region of the parameters of the distributed time-delay controller, and completes the parameterization of the controller. As long as the parameter values are selected in the obtained controller stable set, the multi-delay system can be guaranteed to be stable. Then use the analytical method to obtain the control parameter domain corresponding to the H _∞ performance index, and select parameters in this area to make the system have better dynamic performance and robustness. For the protein synthesis system with time-delay, this method can quickly, effectively and accurately give the control parameter domain of the distributed time-delay controller satisfying the H _∞ performance index, so that the control parameter domain satisfying the H _∞ performance index can be obtained by Select and adjust the parameters to achieve a good control effect.

The present invention is realized through the following technical solutions: use the relevant chemical reactions in the protein synthesis system to carry out system identification on the control object to obtain the model parameters of the controlled object; then calculate the control of the distributed time-delay controller according to the model parameters of the controlled object The boundary line of the parameter stability domain, and judge which side of the boundary line the stability domain is located in, and obtain the stability _domain of the control parameters; The solution method of the stable region calculates the boundary line of the control parameter domain satisfying the H _∞ performance index, and then by judging which side of the boundary line the area satisfying the H _∞ performance index is located, the control parameter domain satisfying the H _∞ performance index is obtained; finally Perform controller tuning, simulation, verification. Selecting control parameters in the obtained control parameter domain satisfying the H _∞ performance index can ensure that the protein synthesis system with time lag has good anti-interference and tracking performance.

Specific steps are as follows:

(1) Before the system enters the design of the distributed time-delay controller, first list the relevant chemical reaction equations of the protein synthesis system, and then use the biological system simulation toolkit in Matlab to convert the relevant chemical reaction equations in the protein synthesis system into normal Differential equations, and then calculate the controlled object model with the following transfer function:

Among them, N(s) and D(s) are polynomials about s, θ represents the time lag, and the highest order of D(s) is greater than the highest order of N(s).

(2) Establish a unit feedback system including a distributed time-delay controller. When the values of the time-delay coefficients τ ₁ , τ ₂ and τ ₃ of the controller are set appropriately, the form of the distributed time-delay controller C(s) is:

Among them, k ₁ , k ₂ and k ₃ are the gain coefficients of each time-delay module in the distributed time-delay controller.

(3) Find the characteristic function of the closed-loop system

(4) For the given τ ₁ , τ ₂ and τ ₃ , by fixing a parameter k ₁ , calculate the boundary line of the stable region composed of the other two control parameters k ₂ and k ₃ on the two-dimensional plane, and obtain the The corresponding stable area is enclosed by the boundary line. The boundary line controlling the parameter space consists of two parts: singular boundary line and non-singular boundary line:

(a) Non-singular boundary line: When s=jω(ω∈(0,+∞)), it is substituted into the closed-loop characteristic function, and the non-singular boundary line is calculated by using the D-segmentation method. The expression of the boundary line is:

in, and Respectively The real and imaginary parts of . Draw the boundary line with values within ω∈(0,+∞).

By judging which side of the boundary line has fewer unstable poles, it is determined which area divided by these boundary lines is the stable domain of the control parameters. In order to determine which side of the non-singular boundary line the stability region of the control parameter is on, the present invention adopts the following method here:

Let s=σ+jω be substituted into the closed-loop characteristic function, and take the real part R and the imaginary part I,

in

because So from the above formula we can get

When k ₁ is fixed, if Q>0 in the two-dimensional plane (k ₂ , k ₃ ), take the direction of the curve along ω∈(0,∞) as the positive direction, and the positive direction of the non-singular boundary line determined by The right side has fewer unstable closed-loop poles than its left side; conversely, the left side of a nonsingular boundary line has fewer unstable closed-loop poles than its right side.

(b) Singular boundary line: when s=0, substitute it into the closed-loop characteristic function to obtain its singular boundary line, the expression of the boundary line is:

k ₃ =-k ₂ -k ₁ (5)

In order to determine which side of the singular boundary line the stable region of the control parameter is located in, the present invention uses the following method to determine:

When s infinitely approaches 0, and Taylor expansion at s=0 has: and respectively substituted into the closed-loop characteristic function to calculate: s+G(0)[k ₁ +k ₂ +k ₃ -s(k ₁ τ ₁ +k ₂ τ ₂ +k ₃ τ ₃ )]＝0, where G(0 ) = G(s) | _s = ₀ . And then there are The necessary condition for the stability of the system is that the poles of the system are all in the left half plane, that is, s<0. When s<0, it can be obtained from the above formula:

(c) Draw the non-singular boundary line and the singular boundary line. The intersection area formed by the side of the non-singular boundary line that does not contain unstable closed-loop poles and the side of the singular boundary line that does not contain unstable closed-loop poles is the stability of the control parameters area.

(5) Calculate the compensation sensitivity function of the closed-loop system

Weighting function Where W _n (s) and W _d (s) are polynomials about s, find the H _∞ performance index that satisfies ||W(s)T(s,k ₁ ,k ₂ ,k ₃ )|| _∞ <γ The control parameter area of , where the H _∞ performance index γ is a real number greater than 0.

given value, the characteristic function satisfying the H _∞ performance index is obtained as:

Substitute s=jω to solve the boundary line satisfying the H _∞ performance index. The expression of the boundary line is:

Different from the stability analysis, here the curve is drawn within ω∈(-∞,+∞) instead of (0,∞). No matter what value k ₁ , k ₂ , k ₃ take, the closed-loop characteristic equation cannot pass through s=0.

For determining which side of the boundary line the parameter domain satisfying the H _∞ performance index is located, the analysis is the same as the analysis for determining which side of the non-singular boundary line the stable domain is. which is

When k ₁ is fixed, if Q>0 in the two-dimensional plane (k ₂ , k ₃ ), there are fewer unstable closed-loop poles on the right side of the positive direction of the determined boundary line than on the left side; otherwise, the left side There are fewer unstable closed-loop poles than to the right of it. It is judged that the boundary line does not contain the intersection area of the unstable closed-loop pole side and the stable domain of the control parameters, and this area is the control parameter domain where the system meets the H _∞ performance index.

(6) Select the control parameters in the control parameter area where the obtained distributed time-delay controller satisfies the H _∞ performance index, and run the system: first, the target gene is transcribed into mRNA, mRNA by adding an invalid gene sequence to the target gene sequence Time lag is generated in each process of translation into protein and protein synthesis, and then the tracking error is accumulated by adding an adder to realize the integral link, thus realizing a distributed time-delay controller in the biological system. Then, activate the target protein gene to make the target gene express and synthesize the target protein; after the synthesized target protein reaches the expected concentration, it starts to inhibit the expression of the target protein gene, thus forming a closed-loop control system. Comparing the synthetic target protein concentration with the initial set concentration produces a tracking error. The tracking error is input to the distributed time-delay controller, and the distributed time-delay controller performs calculations based on the tracking error and outputs control quantities to control the controlled object. The execution process of the distributed time-delay controller can be expressed by the following calculation formula:

where e(t-τ ₁ ), e(t-τ ₂ ) and e(t-τ ₃ ) are the time lags of the tracking error e(t) delaying the transcription of the target gene into mRNA, the translation of mRNA into protein, and the synthesis of protein The values after τ ₁ , τ ₂ and τ ₃ , k ₁ , k ₂ and k ₃ are the control parameters to be obtained, and the e(t-τ ₂ ) and e(t-τ ₃ ) at time 0 to t are calculated by the adder Accumulate and implement integration, and the controller output u(t) acts on the controlled object, so that the controlled object runs in a stable state and meets the H _∞ performance index.

An advantage of the present invention is that it improves the stability and dynamic performance of protein concentration control systems by actively introducing time lag, which is usually considered a disadvantage. And the use of distributed time-delay controllers realizes the occasions where ordinary controllers (such as PID) are difficult to implement. Moreover, the analytical method is used to directly give the control parameter domain calculation method for the distributed time-delay controller of the controlled system to meet the H _∞ performance index. As long as the control parameters are selected in the control parameter domain, the stability and good performance of the closed-loop system can be guaranteed. dynamic performance.

Description of drawings

Fig. 1 is the working flowchart of adopting the method of the present invention.

Figure 2 is a simple model diagram of the system.

Fig. 3 is a closed-loop control structure diagram used in the design method of the distributed time-delay controller in the present invention.

Fig. 4 is a parameter stable region on the (k ₂ , k ₃ ) plane when k ₁ =0, τ ₂ =0.65 min, τ ₃ =1 min in an embodiment of the present invention.

Fig. 5 selects H _∞ performance index γ=0.8, weight function in the embodiment of the present invention , on the (k ₂ , k ₃ ) plane, the parameter domain that satisfies the H _∞ performance index.

Fig. 6 is the unit step response curve DBC1 when control parameters k ₁ =0, k ₂ =3.5, k ₃ =-3 are selected in the stable domain in the embodiment of the present invention, and k ₁ =0, k ₂ are selected outside the stable domain =3.4, the unit step response curve DBC2 when k ₃ =-2.

Fig. 7 is the unit step response curve DBC1 when the control parameters k ₁ =0, k ₂ =3, and k ₃ =-2.5 are selected to meet the H _∞ performance index in the embodiment of the present invention, and the control parameters are selected within the stability range of the control parameters. The unit step response curve DBC2 when the control parameters k ₁ =0, k ₂ =3.3, k ₃ =-2.4 satisfy the H _∞ performance index.

detailed description

The technical solutions of the present invention will be further described below in conjunction with the accompanying drawings and embodiments.

Firstly, the relevant chemical reaction equations in the system are listed, and the biological system simulation toolkit in Matlab converts the relevant chemical reaction equations in the system into ordinary differential equations, and then calculates the controlled object model with the following transfer function, and determines its model parameters; Then, the model parameters and the set distributed time-delay controller parameters τ ₁ , τ ₂ , τ ₃ and k ₁ are implemented in the biological system, that is, by adding an invalid gene to the target protein gene sequence to delay the transcription of the target protein gene into mRNA and mRNA translation into target protein time delay. Based on the D-segmentation method to solve the stable region of the control parameters, the non-singular boundary is calculated by setting s=jω, and the singular boundary is calculated by setting s=0, and further judgment is made on which side of the boundary line the stable region is located, so as to intuitively give The stable domain of k ₂ and k ₃ control parameters surrounded by boundary lines is drawn. Then according to the given H _∞ performance index γ, The sum weight function W(s) solves the boundary line of the H _∞ performance index, further judges which side of the boundary line is the parameter field satisfying the H _∞ performance index, and intersects with the stable area of the control parameters, and the obtained intersection area is It is the parameter field of the control parameters k ₂ and k ₃ satisfying the H _∞ performance index γ.

Example: The control method proposed by the present invention is used in the control system of Hes1 protein synthesis in neural stem cells. When the controlled object is a stable object and the model is unknown, the control parameter domain of the distributed time-delay controller that can make the system stable and satisfy the H _∞ performance index is designed. As long as the control parameters are selected in this parameter domain, the closed-loop system can be guaranteed to have good anti-interference and tracking performance.

The specific implementation steps are given below:

(1) List the relevant chemical reaction equations in the Hes1 protein synthesis system in neural stem cells. The biological system simulation toolkit in Matlab converts the relevant chemical reaction equations in the system into ordinary differential equations, and then calculates the transfer function with the following Control object model:

(2) Time-delay distributed time-delay control parameter k ₁ =0, there is Set the delay time of Hes1 protein gene transcription into mRNA as τ ₂ =0.65min and the delay time of mRNA translation into target protein as τ ₃ =1min, that is, the distributed time-delay controller is:

(3) Find the closed-loop characteristic function of the system as:

(4) When the control parameter k ₁ =0, calculate the non-singular boundary line and the singular boundary line of the stable domain composed of the other two control parameters k ₂ and k ₃ on the two-dimensional plane, and determine the corresponding stable area.

(a) Substituting s=jω(ω∈(0,+∞)) into the characteristic function δ(s), and using the D-segmentation method to calculate the non-singular boundary line, the expression of the boundary line is:

Based on the above expressions, the non-singular boundary of the stable domain can be obtained in the two-dimensional space (k ₂ , k ₃ ), as shown by the black curve in Fig. 4 .

In order to determine which side of the non-singular boundary line the stable domain is on, first substitute s=σ+jω into δ(s), and then calculate When ω∈[0,2], Q>0, that is, in the two-dimensional plane (k ₂ , k ₃ ), the direction along which the curve changes along ω∈(0,∞) is taken as the positive direction. The stable domain of its control parameters is on the right side of the positive direction of the black curve in Figure 4 at the determined non-singular boundary line.

(b) Substituting s=0 into the characteristic function δ(s) to calculate its singular boundary, the expression of the boundary line is

k ₃ =-k ₂ (15)

In order to determine which side of the singularity boundary line the stable region is on, when s is infinitely approaching 0, the Taylor expansion of e ^-0.65s and e ^-s at s=0 is: e ^-0.65s = 1-0.65s, e ^-s = 1-s, and then substitute into δ(s) respectively to calculate Since the necessary condition for system stability is that the poles of the system are all in the left half plane, that is, s<0, then we have

Therefore, it can be seen that the stable region of the control parameters is located in the intersection area above the straight line k ₃ =-k ₂ and below the straight line k ₃ =1.5-0.65k ₂ , that is, the intersection area above the red line and below the blue line in Figure 4 .

(c) Draw the non-singular boundary line and the singular boundary line. The intersection area formed by the side of the non-singular boundary line that does not contain unstable closed-loop poles and the side of the singular boundary line that does not contain unstable closed-loop poles is the area that can make the system stable control parameter field. In the present invention, the specific control parameter stability domain of the closed-loop system is the gray area in Fig. 4 , that is, when k ₁ =0, the values of k ₂ and k ₃ in the gray area can make the closed-loop system stable.

(5) Calculate the compensation sensitivity function of the closed-loop system

Weighting function Take the parameter γ=0.8 of the H _∞ performance index, and obtain the control parameter domain satisfying the H _∞ performance index of ||W(s)T(s,k ₁ ,k ₂ ,k ₃ )|| _∞ <0.8. in,

take again Then solve the characteristic function that satisfies the H∞ performance index as:

Substitute s=jω(ω∈(-∞,+∞)) into V(s) to obtain the boundary line satisfying the H _∞ performance index. The expression of the boundary line is:

Different from the stability analysis, here the curve is drawn when ω takes a value within (-∞,+∞), instead of taking a value at (0,∞). Use ω∈[-2,2] to draw the boundary line of the H _∞ performance index, which is the green curve in Figure 5. The analysis for determining which side of the boundary line the control parameter domain satisfies the H _∞ performance index is the same as the analysis of which side the stable domain is the non-singular boundary line. Let ω take an aggressive value in the interval [-2,2], When ω∈[-2,0) Q>0, and when ω∈(0,2] Q<0, it can be seen that the parameter field that satisfies the H _∞ performance index is located when ω∈(0,2] The right side of the boundary line of the value. The parameter domain obtained by intersecting the area on the right side of the boundary line with the stable domain of the control parameters is the control parameter domain that satisfies the H _∞ performance index, as shown in the purple area in Figure 5.

(6) Select control parameters in the obtained distributed time-delay controller to meet the _H∞ performance index control parameter area, and run the system: first, the target gene is transcribed into mRNA, mRNA by adding an invalid gene sequence to the Hes1 protein gene sequence Time lag is generated in each process of translation into protein and protein synthesis, and then the tracking error is accumulated by adding an adder to realize the integral link, thus realizing a distributed time-delay controller in the biological system. Then, the Hes1 protein gene is activated, so that the Hes1 protein gene is expressed and the Hes1 protein is synthesized. After the synthetic Hes1 protein reaches the expected concentration, it starts to inhibit the expression of the Hes1 protein gene, thus forming a closed-loop control system. Synthetic Hes1 protein concentration compared with the initial set concentration produced tracking error. The tracking error is input to the distributed time-delay controller, and the distributed time-delay controller performs calculation based on the tracking error and outputs the control quantity to control the controlled object. The execution process of the distributed time-delay controller can be expressed by the following calculation formula:

Among them, e(t-0.65) and e(t-1) are the values of tracking error e(t) delaying the time lag of target gene transcription into mRNA and mRNA translation into protein, respectively, and the e from time 0 to t is calculated by the adder The values of (t-0.65) and e(t-1) are accumulated and integrated, and the controller output u(t) acts on the controlled object, so that the controlled object runs in a stable state and meets the H _∞ performance index. Select the control parameters k ₁ =0,k ₂ =3.5,k ₃ =-3 and k ₁ =0,k ₂ =3.4,k ₃ =-2 respectively inside and outside the stability domain of the control parameters, based on the A unit feedback loop verifies whether the system is stable. Given a unit step signal as the input of the system, the system output response curve shown in Figure 6 can be obtained. The corresponding output curve DBC1 of the control parameters selected in the stable domain is stable, and the output curve DBC2 of the selected control parameters outside the stable domain is unstable, and then intuitively verifies the validity of the obtained control parameter stable domain. Then select parameters k ₁ =0, k ₂ =3, k ₃ =-2.5 and k ₁ =0, k ₂ =3.3, k ₃ in the control parameter domain of the H _∞ performance index and outside the control parameter domain of the H _∞ performance index =-2.4, the corresponding unit step response curves are shown in the curves DBC1 and DBC2 in Figure 7. It is verified intuitively that the parameters selected in the control parameter domain of the H _∞ performance index can make the closed-loop system have better dynamic performance.

Claims (1)

1. A method for controlling protein concentration, which is H of a protein synthesis system with time lag_∞The distributed time lag control method comprises the following specific steps:

(1) before the system enters the design of the distributed time lag controller, a relevant chemical reaction equation of the protein synthesis system is listed, then a biological system simulation tool kit in Matlab is utilized to convert the relevant chemical reaction equation in the protein synthesis system into an ordinary differential equation, and a controlled object model with the following transfer function is further calculated:

$G\left(s\right)=\frac{N\left(s\right)}{D\left(s\right)}{e}^{-\mathrm{\&theta;}s}---\left(1\right)$

wherein N(s) and D(s) are polynomials for s, θ represents time lag, and the highest order of D(s) is greater than the highest order of N(s);

(2) establishing a unit feedback system comprising a distributed time lag controller; setting the skew coefficient tau of the controller appropriately₁、τ₂And τ₃The distributed skew controller c(s) is then of the form:

$C\left(s\right)=\frac{1}{s}(k_1{e}^{-{\mathrm{\&tau;}}_{2}s}+k_2{e}^{-{\mathrm{\&tau;}}_{2}s}+k_3{e}^{-{\mathrm{\&tau;}}_{3}s})---\left(2\right)$

wherein k is₁、k₂And k₃The gain coefficient of each time lag module in the distributed time lag controller;

(3) determining a characteristic function of a closed loop system

$\mathrm{\&delta;}\left(s\right)=1+\frac{1}{s}G\left(s\right)(k_1{e}^{-{\mathrm{\&tau;}}_{2}s}+k_2{e}^{-{\mathrm{\&tau;}}_{2}s}+k_3{e}^{-{\mathrm{\&tau;}}_{3}s})---\left(3\right)$

(4) For a given τ₁、τ₂And τ₃By fixing a parameter k₁Calculating two other control parameters k on two-dimensional plane₂And k₃Forming a stable region boundary line, and calculating a corresponding stable region surrounded by the boundary line; the boundary line of the control parameter space consists of a singular boundary line and a non-singular boundary line:

(a) non-singular boundary line: and when s is equal to j omega (omega belongs to (0, and is +/-infinity)), substituting the s into a closed-loop characteristic function, and calculating a nonsingular boundary line by using a D-segmentation method, wherein the boundary line expression is as follows:

$\left\{\begin{array}{c}{k}_2=-\frac{1}{sin\&lsqb;({\mathrm{\&tau;}}_2-{\mathrm{\&tau;}}_3)\mathrm{\&omega;}\&rsqb;}\{\mathrm{Re}\&lsqb;-\frac{j\mathrm{\&omega;}}{G\left(j\mathrm{\&omega;}\right)}-k_1{e}^{-{\mathrm{j\&tau;}}_1\mathrm{\&omega;}}\&rsqb;sin\left({\mathrm{\&tau;}}_3\mathrm{\&omega;}\right)+\mathrm{Im}\&lsqb;-\frac{j\mathrm{\&omega;}}{G\left(j\mathrm{\&omega;}\right)}-k_1{e}^{-{\mathrm{j\&tau;}}_1\mathrm{\&omega;}}\&rsqb;cos\left({\mathrm{\&tau;}}_3\mathrm{\&omega;}\right)\}\\ k_3=\frac{1}{\sin \&lsqb;({\mathrm{\&tau;}}_2-{\mathrm{\&tau;}}_3)\mathrm{\&omega;}\&rsqb;}\{\mathrm{Re}\&lsqb;-\frac{j\mathrm{\&omega;}}{G\left(j\mathrm{\&omega;}\right)}-k_1{e}^{-{\mathrm{j\&tau;}}_1\mathrm{\&omega;}}\&rsqb;sin\left({\mathrm{\&tau;}}_2\mathrm{\&omega;}\right)+\mathrm{Im}\&lsqb;-\frac{j\mathrm{\&omega;}}{G\left(j\mathrm{\&omega;}\right)}-k_1{e}^{-{\mathrm{j\&tau;}}_1\mathrm{\&omega;}}\&rsqb;\cos \left({\mathrm{\&tau;}}_2\mathrm{\&omega;}\right)\}\end{array}\right.---\left(4\right)$

wherein,andrespectively representThe real and imaginary parts of (c), the boundary line is plotted by taking the value in ω ∈ (0, + ∞);

determining which region divided by the boundary lines is a stable region of the control parameter by judging which side of the boundary lines has fewer unstable poles; in order to determine on which side of the non-singular boundary line the stable region of the control parameter is, the invention herein employs the following method:

substituting s-sigma + j omega into the closed-loop characteristic function, taking a real part R and an imaginary part I,

$\left\{\begin{array}{c}R=k_2{e}^{-{\mathrm{\&tau;}}_2\mathrm{\&sigma;}}cos\left({\mathrm{\&tau;}}_2\mathrm{\&omega;}\right)+k_3{e}^{-{\mathrm{\&tau;}}_3\mathrm{\&sigma;}}cos\left({\mathrm{\&tau;}}_3\mathrm{\&omega;}\right)+F_r=0\\ I=-k_2{e}^{-{\mathrm{\&tau;}}_2\mathrm{\&sigma;}}\sin \left({\mathrm{\&tau;}}_2\mathrm{\&omega;}\right)-k_3{e}^{-{\mathrm{\&tau;}}_3\mathrm{\&sigma;}}\sin \left({\mathrm{\&tau;}}_3\mathrm{\&omega;}\right)+F_i=0\end{array}\right.$

wherein

$J_1=\left[\begin{array}{cc}\frac{\&part;k_2}{\&part;\mathrm{\&omega;}}& \frac{\&part;k_2}{\&part;\mathrm{\&sigma;}}\\ \frac{\&part;k_3}{\&part;\mathrm{\&omega;}}& \frac{\&part;k_3}{\&part;\mathrm{\&sigma;}}\end{array}\right]=-{\left[\begin{array}{cc}{e}^{-{\mathrm{\&tau;}}_2\mathrm{\&sigma;}}\cos \left({\mathrm{\&tau;}}_2\mathrm{\&omega;}\right)& e^{-{\mathrm{\&tau;}}_3\mathrm{\&sigma;}}\cos \left({\mathrm{\&tau;}}_3\mathrm{\&omega;}\right)\\ -e^{-{\mathrm{\&tau;}}_2\mathrm{\&sigma;}}\sin \left({\mathrm{\&tau;}}_2\mathrm{\&omega;}\right)& -e^{-{\mathrm{\&tau;}}_3\mathrm{\&sigma;}}\sin \left({\mathrm{\&tau;}}_3\mathrm{\&omega;}\right)\end{array}\right]}^{-1}\left[\begin{array}{cc}\frac{\&part;R}{\&part;\mathrm{\&omega;}}& \frac{\&part;R}{\&part;\mathrm{\&sigma;}}\\ \frac{\&part;I}{\&part;\mathrm{\&omega;}}& \frac{\&part;I}{\&part;\mathrm{\&sigma;}}\end{array}\right]$

Due to the fact thatSo that can be obtained from the above formula

$Q=Sgn\left(\frac{\&part;k_2}{\&part;\mathrm{\&omega;}}\frac{\&part;k_3}{\&part;\mathrm{\&sigma;}}-\frac{\&part;k_2}{\&part;\mathrm{\&omega;}}\frac{\&part;k_3}{\&part;\mathrm{\&sigma;}}\right)=Sgn\{sin\&lsqb;({\mathrm{\&tau;}}_3-{\mathrm{\&tau;}}_2)\mathrm{\&omega;}\&rsqb;\}$

When k is₁When fixed, if Q > 0 two-dimensional plane (k)₂,k₃) Taking the direction of the curve changing along omega ∈ (0, ∞) as the positive direction, and determining that the positive right side of the nonsingular boundary line has fewer unstable closed-loop poles than the left side of the nonsingular boundary line;

(b) singular boundary line: when s is equal to 0, substituting the s into a closed-loop characteristic function to obtain a singular boundary line of the closed-loop characteristic function, wherein the boundary line expression is as follows:

k₃＝-k₂-k₁(5)

in order to determine on which side of the singular boundary line the stable region of the control parameter is located, the present invention makes a decision by using the following method:

as s approaches to 0 at infinity,andtaylor expansion at s-0 has: and substituting the closed-loop characteristic functions respectively to obtain: s + G (0) [ k ]₁+k₂+k₃-s(k₁τ₁+k₂τ₂+k₃τ₃)](ii) 0, wherein G (0) ═ G(s) does not smoke_s＝0(ii) a And then haveThe necessary condition for system stability is that the poles of the system are all in the left half-plane, i.e. s < 0, when s < 0, we can see the above equation:

$\left\{\begin{array}{c}{k}_1+k_2+k_3>0\\ k_1{\mathrm{\&tau;}}_1+k_2{\mathrm{\&tau;}}_2+k_3{\mathrm{\&tau;}}_3<\frac{1}{G\left(0\right)}\end{array}\right.---\left(6\right)$

(c) drawing a non-singular boundary line and a singular boundary line, wherein an intersection region formed by one side of the non-singular boundary line, which does not contain the unstable closed-loop pole, and one side of the singular boundary line, which does not contain the unstable closed-loop pole, is a stable region of the control parameters;

(5) calculating a compensated sensitivity function for a closed loop system

$T(s,k_1,k_2,k_3)=\frac{N\left(s\right)(k_1{e}^{-{\mathrm{\&tau;}}_{1}s}+k_2{e}^{-{\mathrm{\&tau;}}_{2}s}+k_3{e}^{-{\mathrm{\&tau;}}_{3}s})e^{-\mathrm{\&theta;}s}}{sD\left(s\right)+N\left(s\right)(k_1{e}^{-{\mathrm{\&tau;}}_{1}s}+k_2{e}^{-{\mathrm{\&tau;}}_{2}s}+k_3{e}^{-{\mathrm{\&tau;}}_{3}s})e^{-\mathrm{\&theta;}s}}$

Weighting functionWherein W_n(s) and W_d(s) is a polynomial expression for s, and T (s, k) satisfying | | W(s)₁,k₂,k₃)||_∞< gamma H_∞Control parameter region of performance index, where H_∞The performance index γ is a real number greater than 0;

$W\left(s\right)T(s,k_1,k_2,k_3)=\frac{W_n\left(s\right)}{W_d\left(s\right)}\frac{N\left(s\right)(k_1{e}^{-{\mathrm{\&tau;}}_{1}s}+k_2{e}^{-{\mathrm{\&tau;}}_{2}s}+k_3{e}^{-{\mathrm{\&tau;}}_{3}s})e^{-\mathrm{\&theta;}s}}{sD\left(s\right)+N\left(s\right)(k_1{e}^{-{\mathrm{\&tau;}}_{1}s}+k_2{e}^{-{\mathrm{\&tau;}}_{2}s}+k_3{e}^{-{\mathrm{\&tau;}}_{3}s})e^{-\mathrm{\&theta;}s}}---\left(7\right)$

given aThe value of (A) is solved to satisfy H_∞The characteristic function of the performance index is:

substituting s as j omega to solve the problem of H_∞A boundary line of the performance indicator, the boundary line having an expression:

unlike the stability analysis, the curve is plotted by taking the value in ω ∈ (-infinity, + ∞) instead of taking the value in (0, ∞), regardless of k₁,k₂,k₃Taking the value, the closed-loop characteristic equation cannot pass through s being 0;

satisfies H for the determination_∞The parameter domain of the performance index is positioned on which side of the boundary line, and the analysis of the parameter domain is the same as that of which side of the non-singular boundary line of the stable domain is judged; namely, it is

$Q=Sgn\left(\frac{\&part;k_2}{\&part;\mathrm{\&omega;}}\frac{\&part;k_3}{\&part;\mathrm{\&sigma;}}-\frac{\&part;k_2}{\&part;\mathrm{\&omega;}}\frac{\&part;k_3}{\&part;\mathrm{\&sigma;}}\right)=Sgn\{sin\&lsqb;({\mathrm{\&tau;}}_3-{\mathrm{\&tau;}}_2)\mathrm{\&omega;}\&rsqb;\}$

When k is₁When fixed, if Q > 0 two-dimensional plane (k)₂,k₃) The positive right side of the determined boundary line has fewer unstable closed-loop poles than the left side, and the left side has fewer unstable closed-loop poles than the right side; judging the intersection region formed by the side of the boundary line not containing the unstable closed loop pole and the stable region of the control parameter, wherein the intersection region is the region which satisfies the requirement H of the system_∞A control parameter field of the performance indicator;

(6) satisfying H in the obtained distributed time lag controller_∞Selecting control parameters from the control parameter area of the performance index,operating the system: firstly, adding an invalid gene sequence into a target gene sequence to enable the target gene to be transcribed into mRNA, translate the mRNA into protein and generate time lag in each process of synthesizing the protein, and adding an adder to enable tracking errors to be accumulated to realize an integration link, thereby realizing a distributed time lag controller in a biological system; then, activating the target protein gene to express the target gene and synthesize the target protein; after the synthesized target protein reaches the expected concentration, the expression of the target protein gene is inhibited, so that a closed loop control system is formed; comparing the concentration of the synthesized target protein with the initial set concentration to generate a tracking error; the tracking error is input to a distributed time-lag controller, the distributed time-lag controller carries out operation based on the tracking error and outputs a control quantity to control a controlled object; the execution process of the distributed skew controller can be represented by the following calculation formula:

$u\left(t\right)=k_1{\&Integral;}_0^{t}e(t-{\mathrm{\&tau;}}_1)dt+k_2{\&Integral;}_0^{t}e(t-{\mathrm{\&tau;}}_2)dt+k_3{\&Integral;}_0^{t}e(t-{\mathrm{\&tau;}}_3)dt---\left(10\right)$

wherein e (t- τ)₁)、e(t-τ₂) And e (t- τ)₃) Delay lag tau through transcription of target gene into mRNA, translation of mRNA into protein and protein synthesis, respectively₁、τ₂And τ₃Value of latter, k₁、k₂And k₃For the control parameter to be found, e (t-tau) from 0 to t is added by an adder₂) And e (t- τ)₃) Accumulating to implement integration, and applying the output u (t) of controller to the controlled object to make it operate in stable state and meet H_∞Performance index.