1. Introduction
Over the recent decades, there had been increasing interest in control theoretic approach in synthetic biology [1]. Especially when dealing biomolecular circuits, external disturbances caused from molecular events are present which impose difficulties in observation of the state variables and design of the feedback controllers [2]. To attain certain amount of robustness and performance of the feedback control, quantitative inspection of the synthetic gene circuits was attempted. Regarding the traditional proportional-integral-differential (PID) control design, [3] demonstrated for gene circuits under stochastic disturbance. The biggest merit of PID control design is that the priori knowledge about the system is not necessary for the design of feedback control. In variational approaches, [4] used bang-bang control, which is typically obtained by solving the maximum principle. As there were animated researches held with the model predictive control (MPC) recently, [5] compared the three approaches.
Meanwhile at the other side of the academia, there had been breakthroughs on numerical solution methodology of nonlinear partial differential equations (PDEs), namely the Galerkin’s Approximation. The approximation is done by finding an appropriate linear combination of generating functions, or basis functions, to reduce the size of the residual of PDE to zero. Such method is meaningful since the typical optimal control problems are involved with the nonlinear PDEs, the Hamilton-Jacobi-Bellman (HJB) and Hamilton-Jacobi-Isaacs (HJI) equations. By solving these PDEs, the optimal control input required for minimization of the performance measure is determined. The PDE that needs to be solved depends on the performance measure. For typical choices so-called the
and
feedback control, HJB and HJI need to be solved respectively [6]. But when the system is nonlinear, these PDEs are well-known to be extremely hard to solve. By exploiting Galerkin’s approximation, [7] had established a rigorous mathematical approach by using Galerkin’s approximation in iterative manner to solve HJB. Such methodology was named Successive Galerkin Approximation (SGA). Then it was proposed in [8] a similar approach to solve HJI and the point-wise convergence properties were also proved. Performance of such method was also verified. Moreover, it was confirmed in [9] about the superior performance of the controller designed through SGA over controller designed by linearization of the nonlinear dynamics. Although such methodology also requires precise modeling of the system we aim to control, thanks to chemists and biologists, a large portion of biochemical systems have reliable state dynamics and stochastic modeling of the intrinsic disturbances [10]. Unlike the PID control,
-feedback control provides quantitative amount of margin of the control system that can withstand against the external disturbance which makes it possible to maximize the disturbance rejection performance. Hence such optimal control strategy is preferred over traditional control strategies if an accurate modeling of the system is available.
To this end, we present an implementation of
and
feedback control design methodology with SGA. Along the simulation, several limitations of the SGA algorithms are indicated and mathematical reasonings are provided.
2. Feedback Control Design Problem
Let us consider the smooth nonlinear control system represented in the below form.
(1)
where
,
,
for each t, and
,
,
.
denotes the state feedback. The system is observable through output
where
is a smooth function with
. When the unknown external disturbance
is entering the system, it is desired to find an appropriate state feedback
that achieves robust stability and performance. Throughout this paper, we impose the below assumption about the system. With this assumption, we introduce the two design purposes, namely
and
-feedback throughout this section.
Assumption 1
and
.
2.1.
-Feedback Control Design
For
-feedback control design problem, the performance measure is set as
where
,
for some
. It was shown from the previous studies including [11] that for the solution
of the Hamilton-Jacobi-Bellman Equation (HJB) below, feedback control
attains minimum of the above performance measure.
(2)
Although the PDE (2) can be easily solved by casting into form of algebraic Riccati equation in linear cases, gene networks typically include high nonlinearities and hence PDE (2) is extremely hard to solve in general.
2.2.
-Feedback Control Design
Unlike the
-feedback control, the objective of
-feedback control is to achieve asymptotic stability and disturbance-to-output
-gain under certain level
. Such condition can be expressed as below.
(3)
where
and
for some
. It was proved in [6] that for solution
of the Hamilton-Jacobi-Isaacs Equation (HJI) below, feedback control
achieves the condition of Equation (3). Moreover, it was also proven that the worst case disturbance if given by
. The inequality of Equation (3) is tight for such choice of
and
and when
.
(4)
Like the HJB, HJI is not tractable in general cases.
3. Successive Galerkin Approximation
A typical strategy to solve PDEs in numerical manner is Galerkin method [12]. Galerkin method approximates the solution with linear combination of finite number of trial functions.
[8] applied Galerkin method in sequential manner to obtain convergent solutions of HJB and HJI. Point-wise convergence of such method was also proved. Such approach was named Successive Galerkin Approximation (SGA). Although convergence and consistency under finite number of trial functions was not provided, with reasonable number and choice of trial (or basis) functions, it was verified that the feedback design through SGA performs decently. By taking the notations from [8], we summarize the SGA algorithm for HJB and HJI respectively in this section.
Let us begin by introducing the notations. For the basis functions
, let
. Now when the solution is approximated by
, let
. We further denote the integration region of Galerkin method by
. Now we begin by introducing the iteration matrices for solution of the HJB.
(5)
The derivation of the above formulas (5) can be found at [8].
In Algorithm 1,
is the result from the
-th iteration and
is the tolerance for termination criterion.
For solution of HJI, construction of the iteration matrices is as below.
(6)
In Algorithm 2,
is the result from the
-th iteration, or the i-th inner loop and j-th outer loop iteration.
is the tolerance for termination criterion. It was proved in the previous researches the convergence of both Algorithm 3 and Algorithm 3 for complete set of basis functions of continuous
space. But under finite number of basis functions, convergence of the solution procedure and convergence of the residual of the nonlinear PDE to zero was not proved. Due to finite computation resources, this is the case in practical applications. Our consequent results in this paper will review the results obtained from the SGA algorithms and point out such issues especially for the application in HJI case.
4. Problem Description
For illustrative purpose, we consider the gene regulatory network in cascaded form brought from [10]. The system dynamics
is represented by the below first order ODE. Though we specify the system dynamics as below, large amount of gene network systems are represented in the similar form and hence such an approach can be employed at various situations [10]. For the
-feedback control design of this system, linearization and application of fuzzy interpolation was attempted in [2].
(7)
where
is the external disturbance and the initial conditions are given by
,
,
,
. Such representation is also called the generalized-mass-action(GMA) form. It is evident that one equilibrium point is
. We sketch a formulaic design method of
and
-feedback controls for this system. Rewriting the Equation (2), it follows
(8)
Because boundness of
has to be guaranteed for Galerkin’s approximation to be used,
. Hence it is evident that it is required to have
for application of the SGA algorithm. Due to such reason, we consider the perturbed system as below.
(9)
By such perturbation, the system now has
as an equilibrium point.
5. Numerical Demonstration & Discussion
In this section, we provide the numerical demonstration results of the previously described nonlinear gene network system. For both
and
-feedback designs, total six basis functions were used and are given by
. Such choice was to obtain linear approximation of the true solution
of HJB or HJI and satisfy
simultaneously. The integration region was set as
. For the feedback input,
was considered where
and
are the 2 × 2 identity and zero matrices respectively. Numerical simulation was done by Simulink software on Matlab 2019b with Macbook Pro, 2.4GHz Intel Core i5.
5.1.
-Feedback Design
For the design parameters,
and
was chosen. With the given basis, it took 7 iterations to converge within the tolerance
. Converged solution is shown in Table 1. The state and control variables from the designed
controller are depicted in Figure 1. Simulation results with external disturbance entering through the state dynamics of
are depicted in Figure 2 together with the random disturbance. The designed control system withstands such disturbance decently.
(a) State variables(b) Control variables
Figure 1. Time profiles of the state variables with the designed
-feedback: without external disturbance.
(a) State variables(b) Control variables(c) Applied External Disturbance
Figure 2. Time profiles of the state variables with the designed
-feedback: with external disturbance.
Table 1. Obtained solution of the HJB.
5.2.
-Feedback Design
For the
-feedback control design, as from [2], it is desirable to find the smallest
that the solution of HJI exists. Hence for the design parameter,
was used for the initial choice and we designed the feedback control throughout the solution. It took total 8 iterations for the outer loop to converge while inner loops typically took 4 iterations to converge within the same tolerance from previous subsection. Converged solution is shown in Table 2. Though the
feedback controller had shown slightly rapid convergence to the desired equilibrium point, both controllers show similar performances. Comparing with the results from [2] obtained from fuzzy interpolation after linearization, Figure 3 and Figure 4 show more rapid convergence. Such superiority of the feedback design without linearization was also confirmed at [9]. Then by trial and error, it was confirmed that the SGA algorithm also converges for smaller values of
. But for the solution obtained for
less than 2, the SGA algorithm did converge but the designed feedback controller exhibited unstable properties. Converged solution is shown in Table 3. Simulation results for
is depicted in Figure 5 and Figure 6 which shows unstable response regardless the presence of disturbance. We speculate that this is due to non-zero residual of the HJI (4). Because the choice of our basis was obviously not a complete set that can span the whole continuous
space, although the projection of the residual onto the subspace generated by our basis functions may be zero, it does not imply that the residual itself is zero. Regarding this, we would like to point out that the solution obtained from the SGA algorithm is indeed the utmost solution with the given basis, but it does not necessarily imply that it is the exact solution of HJI. This could be resolved through adding the number of basis functions or choosing appropriate basis functions depending on the system dynamics through heuristic. Though it was proposed in [8] a clever methodology exploiting tensor products to mitigate the computational burden of massive number of numerical integrations, in practice we are only allowed to use finite number of basis functions due to finite computing power. Also it should be pointed out that in a lot of cases, full-state observability of the system is not present and hence there are limitations in choice of basis functions too.
Table 2. Obtained solution of the HJI with
.
(a) State variables(b) Control variables
Figure 3. Time profiles of the state variables with the designed
-feedback (
): without external disturbance.
(a) State variables(b) Control variables(c) Applied External Disturbance
Figure 4. Time profiles of the state variables with the designed
-feedback (
): with external disturbance.
(a) State variables(b) Control variables
Figure 5. Time profiles of the state variables with the designed
-feedback (
): without external disturbance.
Table 3. Obtained solution of the HJI with
.
(a) State variables(b) Control variables(c) Applied External Disturbance
Figure 6. Time profiles of the state variables with the designed
-feedback (
): with external disturbance.
6. Conclusions & Future Scope
In this paper, nonlinear
and
-feedback controls have been developed for nonlinear gene regulatory network in GMA form. With the choice of basis functions from Section 5, Successive Galerkin Approximation did converge rapidly for both the Hamilton-Jacobi-Bellman and Hamilton-Jacobi-Isaacs equations. Because the majority of gene network systems are of the form of Equation (7), such feedback design methodology from this paper could be applied in a wide range of systems in GMA form. Comparison between the designed
and
-feedback control systems was made and for the hyper parameters used, it was confirmed that the two controllers exhibited similar performances.
It was also confirmed from our simulation that even when the SGA algorithm converged, there were cases where the corresponding feedback control system was unstable. This is due to nonzero residual of the HJI, and it needs to be further resolved by improvements in the algorithm for solving HJB and HJI. This issue could indeed be overcome by different choice of basis functions, but still there is no guarantee that the solution obtained with new basis functions will not exhibit such instability. Additionally in practice, there are limitations in observability of the state variables, which again impose restriction on the choice of the basis functions. Hence it requires further research to analyze the feedback control systems’ properties under finite number of basis functions to overcome these issues. Our consequent work will include breakthroughs in such issues discussed and consideration of stochastic disturbance entering through the system dynamics due to molecular noises.