Computational optimal control of 1D colloid transport by solute gradients in dead-end micro-channels

Diffusiophoresis is a common phenomenon that occurs when colloids are placed in the non-uniform solute concentration. It generates solute gradients which force the colloids to transfer toward or away from the higher solute concentration side. In this paper, we consider the input sequence control of the colloid transport in a dead-end micro-channel with a boundary solute concentration being manipulated, which has a wide range of applications such as drug delivery, biology transport, oil recovery system and so on. We model this process by a coupled system, which involves the solute diffusion equation and the colloid transport model. Then an optimal control problem is formulated, in which the goal is to minimize colloid density distribution deviation between the computational one and the target at a pre-specified terminal time. To solve this partial differential equation (PDE) optimal control problem, we first apply the control parameterization method to discretize the boundary control and transfer it into an optimal parameter selection problem. Then, using the variational method, the gradient of the objective function with respect to the decision parameters can be derived, which depends on the solution of the coupled system and the costate system. Based on this, we propose an effective computational method and a gradient-based optimization algorithm to solve the optimal control problem numerically. Finally, we give the simulation results to demonstrate that the objective function based on the proposed method is less nearly two orders of magnitude than that of a constant value control strategy, which well illustrates the effectiveness of the proposed method.


1.
Introduction. Colloid transport propelled by solute gradients has attracted much scientific interest recently [44,20]. This phenomenon, also called diffusiophoresis [3], commonly occurs in the solute-particle system such as drug delivery [35], biology transport [2], oil recovery system [10] and so on. The size of colloid usually ranges from 10nm to 10µm, less than the size of pipe channels [48,17]. The physical mechanism mainly originates from the interaction between colloid's surface and the solute. The non-uniform solute concentration generates solute gradients, which lead to the active motion of colloid [19,18]. Under such circumstances, the transport rate of colloid is much higher than pure diffusion.
Dead-end channels are commonly found in many areas. Colloid transport into such dead-end geometries efficiently is important. In drug delivery process, transferring the drug into blood dead-end pores is important for the recovery of the patient, which greatly improves the effectiveness of the drug absorption [35]. In biological transport, organisms such as bacteria or cells move towards the nutrients or poisons in the dead-end, which promotes their decomposition or breakdown [2]. In oil refinery, pumping the colloidal particles and oil emulsions in and out of the dead-end capillaries is significant implications for oil recovery systems [19,18].
Generally, the way of colloid transport can be divided into two categories: external induction and interior induction. For the external induction, colloid transport is driven by external force such as electric fields [29], magnetic fields [18], thermal fields [42] and so on. The induction depends on the harsh external conditions and sometimes limits the applications on a large scale. While in the interior induction, colloid transport is usually driven by interior inductions such as chemical reaction [21], osmophoresis [38] and diffusiophoresis [41]. The induction is less dependent on external conditions and a good supplement for the external induction. In practice, we often combine two driven categories to generate better results. In this paper, we pay attention to one of the interior inductions: solute gradient.
The mathematical equations describing the solute concentration distribution along the long straight pipe are well studied. Callewaert [12] considers the numerical constant of the Taylor-Aris expression for the dispersion in a laminar flow and uses the MATLAB PDE toolbox for calculation. Wu [50] considers the longitudinal normality of mean concentration and explores the approach towards transverse uniformity of the solute concentration distribution. Kovalchuk [5] deals with the description of stationary transport of non-electrolyte molecules or electro-neutral particles due to hydrodynamic dispersion. Abderrezzak [1] considers one-dimensional (1-D) numerical models of solute transport with unknown parameters in streams under steady and unsteady flow conditions. In this paper, since the diffusion process of solute concentration is dominant, we consider a reduced diffusion equation. As for the mathematical equations describing the colloid transport model, we apply the classical advection-diffusion equation proposed in [25].
The diffusiophoresis process, which involves the solute diffusion equation and colloid transport model, provides physical mechanisms for advanced control strategy study with physical manipulated boundary actuation, such as the solute concentration and colloid density. However, few work of advanced control strategies for this process are studied [6]. In this paper, we propose a PDE-constrained optimization problem applied to the colloid transport with a solute concentration control on the boundary, which greatly improves the flexibility of colloid transport. But the direct design of optimal control strategies for the coupled PDEs is usually challenging since analytical solutions are hard to obtain if it is not impossible. Thus, effective computational optimal algorithms are necessary to solve the optimal control problem of coupled PDEs.
Optimization and controller synthesis based on PDE models have a long tradition in many applications. Also, with the advances of sensor and actuation technologies, process control via PDE models becomes feasible and there are quite a few successful stories in applications, including soaking pit [37], thermo-acoustic combustion instability control [4], aerodynamic flows [40], fuel cells [8], etc., just to name a few.
The paper is organized as follows. In Section 2, we introduce a coupled PDE system to describe the colloid transport in the dead-end micro-channel. Then, we formulate a PDE optimal control problem for the colloid transport in the dead-end micro-channel (as shown in Figure 1), where the inlet solute concentration is manipulated. In Section 3, we first discretize the inlet boundary solute concentration by using the control parameterization method and convert it into an optimal parameter selection problem. Then, we use the variational method to derive the gradient of the objective function, which depends on the solution of the coupled system and the costate system. In Section 4, we propose an effective computational method and a gradient-based optimization algorithm to solve the optimal control problem arising in colloid transport. In Section 5, numerical simulation is given to show effectiveness of the proposed computational optimal control strategy by comparing the target colloid density distribution and the computational one.
2.1. Mathematical model of the colloid transport. We now consider the situation where the entrance of a dead-end micro-channel of the length L = 400µm is connected to a main channel as shown in Figure 1 and consider the colloid transfer in the dead-end micro-channel [46]. The main channel, which is full of NaCl solution carrying the colloids, transports solute from end to end at a certain speed. An oil droplet or an air bubble separates well-mixed solute into different chambers in the main channel. These designs make a given admissible boundary control sequence possible for the dead-end micro-channel [27,16]. In practical manipulation, we can design the micro-channel structure combined with the droplet microfluidics technique, which can generate an oil droplet or an air bubble. The micropump can transport the microfluid along the main channel. Compared to the main channel, the size of dead-end micro-channel is far less. Thus, well-mixed solute in each separated chamber and colloid density connected to the dead-end micro-channel can be seen as a constant. By neglecting the effects of temperature variation and solute separators, the solute diffusion equation in the dead-end micro-channel can be described by following dimensionless diffusion equation [50], where C(x, τ ) is the dimensionless solute concentration, x ∈ [0, 1] the spatial variable along the dead-end micro-channel, τ = Dst L 2 the dimensionless time, t ∈ [0, T ], T the total time span and D s the ambipolar diffusion coefficient for solute. The boundary conditions at the entrance of the dead-end micro-channel are piecewise constant and the dead-end of it satisfies a natural boundary condition. Thus, the where u(τ ) denotes the boundary inlet solute concentration, C(1, τ ) = 0 the solution at x = 1 always exists. The boundary solute concentration u(τ ) in the dead-end channel has the direct physical meaning and can be regulated combined with the experimental setup as shown in Figure 1. The solute gradient can not be regulated directly by the microfluidic devices and should be regulated through the solute concentration in practice. Moreover, the solute concentration can be easier measured by the planar laser-induced fluorescence technique while the solute gradient can not be measured directly. Thus, the solute concentration condition chosen as the boundary control variable is more reasonable. Furthermore, the initial condition for system (1) is where Ψ(x) is a given function describing the initial state of the solute concentration along the dead-end micro-channel. The transient nature of the solute diffusion equation generates solute gradients, which drive the colloids into or outward the dead-end micro-channel from the main channel. We assume that the solute gradients can be linearized on the scale of the colloids and obtain the colloid velocity along the axial velocity direction as follows [3] V Thus, the advection-diffusion equation of coupled dimensionless colloid transport model is where N (x, τ ) is the dimensionless colloid density, D c the colloid diffusivity, Γ c the diffusiophoretic mobility of a colloid. The first term of the right-hand-side in (5) denotes the colloid diffusion transport while the second term in (5) denotes the colloid transport by the propellant of the solute gradient, which is dominant during the colloid transport process. For more information on the colloid transport model, please refer to the supporting information of [46] and [25]. The boundary conditions of system (5) are whereN is a constant colloid concentration. The initial condition is where Φ(x) is a given function describing the initial state of the colloid density along the dead-end micro-channel.
2.2. The optimal boundary control problem. Our goal is to choose an admissible boundary control u(τ ) so that the actual colloid density distribution can be brought as close as possible to the target colloid density distribution at the terminal time (τ = Z, where Z = DsT L 2 ). Thus, the objective function can be given as follows: where N d (x) is the target colloid density distribution along the space at the terminal time. The objective function (8) penalizes the deviation between the actual colloid density distribution and the target colloid density distribution. Its aim is to transport the colloids including the peak position into the dead-end efficiently and reasonably at a pre-specified time.
The inlet boundary solute concentration u(τ ) is required to satisfy the following bound constraint: where u max denotes the maximum control parameter. Any control u(τ ) that satisfies the bound constraint (9) is called an admissible boundary control policy. Let U be the class of all such admissible control. Our optimal boundary control problem is now defined as follows.
Problem P 0 . Given the mathematical model (1)-(7), we choose the boundary control u(τ ) ∈ U to minimize the objective function (8) subject to the bound constraint (9).
Problem P 0 is a nonlinear optimal control problem of infinite dimensional systems governed by PDEs. Recently, there have been mainly three solutions for nonlinear optimization problem, including control vector parameterization, spatial discretization and orthogonal collocation. For the control vector parameterization, the finite difference approach, state trajectory sensitivity approach and adjoint system approach for ordinary differential equations (ODEs) to obtain the gradients are examined [43]. The control parameterization method for nonlinear optimal control of ODEs is investigated [30,51]. For the spatial discretization, the method of lines [45], which reduces the PDEs to ODEs, is proposed.
For the orthogonal collocation, an effective orthogonal collocation method is developed by introducing Lagrange interpolation polynomial. Thus, the differential equation model is reduced to an equality constrained nonlinear programming (NLP), which can be solved by the specialized methods [7]. However, this process leads to a larger nonlinear programing problem.
Our novelty is summarized as follows. Firstly, we focus on the control-oriented modelling and optimal control of the colloid transport in microfluidic systems. This work is seldom studied but very promising, which promotes the combination of informatics and microfluidics. Secondly, the solutions of optimal control problem mentioned above ( [7,28]) belong to the "discretize-then-optimize" category, which reduces the PDEs to ODEs or algebraic equations. Thus, many optimal control methods can be applied. In this paper, we keep the PDE structure of the colloid transport system but only discretize the control input function. For the numerical simulation of the original system and costate system, we can use various effective numerical calculation methods, even computational packages (e.g., COMSOL Multiphysics [22], FINICS [34], etc.).
3.1. Control parameterization. To solve Problem P 0 , we first subdivide the time donates the kth control subinterval and τ 0 = 0, τ p = Z [47]. Since the inlet solute concentration u(τ ) is piecewise constant, we can discretize the control as where σ k is the value of u(τ ) on each control subintervals [τ k−1 , τ k ). We can rewrite (10) as where We denote σ = σ 1 , . . . , σ p ∈ R p , then the objective function (8) becomes where N p (x, τ ) denotes the solution of system (1)- (7) with u(τ ) = u p (τ ), so as C p (x, τ ). Furthermore, the bound constraint (9) becomes We define and let U p be the set of all piecewise constant function given in the form of (10) with σ k , k = 1, . . . , p satisfying the constraints (14). Clearly, each u(τ ) ∈ U p corresponds to a vector σ ∈ U uniquely. Thus, Our optimal parameter selection problem is stated as follows.

Gradient computation.
Problem P 1 is a PDE-constrained nonlinear optimization problem. In order to solve it from the view of computational optimal control, we need the gradient of the objective function (13) with respect to the control parameter vector σ [49]. However, the objective function (13) depends on σ implicitly through the coupled system (1)- (7), computing its gradient is a nontrivial task. In this paper, we use the variational methods to obtain the formulas of the required gradient.
Theorem 3.1. The gradient of the objective function in Problem P 1 is given by where α(0, τ ), β(0, τ ) are obtained by solving the following costate system with the boundary conditions, and the terminal conditions Concrete proof is given in Appendix A.
Remark 1. The colloid transport model contains the strong nonlinear term (the second term of the right-hand-side in (5)). Due to the nonlinear term ln(C), the optimal control problem defined in Problem P 1 is lack of the Frechet differentiability in L p . Thus, the derivation of the gradient in this paper is only a formal calculation.

4.
Numerical solution procedure. Based on the gradient derived in Section 3, we can propose the following conceptual framework for solving Problem P 1 using standard optimization packages, such as FMINCON.
Note that Steps 6-9 can be implemented automatically by existing nonlinear optimization solvers, such as FMINCON in MATLAB. Since we have the PDE model of the coupled system and the costate system, effective numerical computating is still needed. For dimensionless diffusion equation (1), we use implicit difference scheme to approximate and obtain where 0 ≤ ϑ ≤ 1, O(∆τ ) and O(∆x 2 ) denote omitted first and second order terms such that O(∆τ ) → 0 as ∆τ → 0 and O(∆x 2 ) ∆x → 0 as ∆x → 0. From the boundary conditions (2), we have From the initial condition (3), we have Similarly, we obtain the following approximation for the dimensionless colloid transport model (5). As for ∂N /∂x, we use the upwind scheme [24], From the boundary conditions (6), we have From the initial condition (7), we have

4.2.
Numerical implementation of the costate system. We also use the finitedifference method to solve the costate system (17) and (18). Using Taylor expansion and the upwind scheme, we obtain From the boundary conditions (19), we have From the initial condition (20), we have Furthermore, for the numerical integration of the objective function (8) and the gradient (16), we apply the composite Simpson's rule. More details on numerical integration algorithms are available in [11].
From the simulation results, the spatiotemporal evolution of colloid transport by the finite difference method well approximates the results in [46] with the same physical parameters and conditions. It shows the stiff problem for the coupled PDEs of consideration is not serious and errors by using the finite difference method are acceptable in this case. Also, to improve the computational accuracy, the computational step size can be taken to be sufficiently small. However, this leads to the longer computational time. In order to make the computation more efficient, we can apply some numerical methods to deal with the stiff problem. For example, the method of lines [45] can be utilized to reduce the PDEs to ODEs. Then, implicit-explicit methods, such as the implicit-explicit linear multistep methods and the implicit-explicit Runge-Kutta methods [26] can be applied to handle the computational accuracy issue of stiff equations. 5. Numerical simulations. For numerical simulations, we set parameters of the solute diffusion equation: ambipolar diffusion coefficient for solute D s = 1.61 × 10 3 µm 2 /s, initial state of the solute concentration Ψ(x) = 2mmol/L and set parameters of the colloid transport model: diffusiophoretic mobility of a colloid Γ c = 610µm 2 /s, colloid diffusivity D c = 0.48µm 2 /s, boundary conditionN = 1 × 10 −3 µm −3 and Φ(x) = 0 (there is no colloids in dead-end micro-channel initially) for the colloid diameter of 1.01µm, which are consistent with those in [46]. We also set the number of space and time subintervals as M X = 50, M T = 6000 in the numerical computation, respectively.
Since the gradient ∂ ln C(x,τ ) ∂x usually decreases along the dead-end micro-channel, the colloid velocity decreases as the colloids move deeper into the micro-channel accordingly. Thus, colloids accumulate near the leading edge of the colloid transport, which leads to the "colloid focusing" effect. By assuming the parts to be treated ranges from 70% to 85% of the dead-end channel, we choose the target colloid density distribution N d (x) involved in (8) as shown in Figure 2(a) where the peak position of the colloid is around 80%. The upper bound constraint in (9) is chosen as u max = 2, which is equal to the initial solute concentration in the dead-end micro-channel. Under the parameters of microfluidic system above, the peak position of the colloids is around 80% within 100s − 300s, so we choose the total time span as T = 200s reasonably.
The numerical simulation was carried out within the MATLAB programming environment (version R2010b) running on a personal computer with the following configuration: Intel Core i5-2320 3.00GHz CPU, 4.00GB RAM, 64-bit Windows 7 Operating System. We apply our MATLAB code to implement the gradient-based optimization procedure as shown in Section 4 Algorithm 1 by using FMINCON (to perform the optimization steps).
We subdivide the time horizon [0, Z] into p = 4 control subintervals [τ k−1 , τ k ), k = 1, . . . , 4 with τ 0 = 0 and τ 4 = Z. Thus, the control input u(t) is parameterized as σ k , k = 1, ..., 4. The initial guess of parameterized boundary control is chosen as σ k = 1, k = 1, ..., 4. Figure 2(b) shows the optimal parameterized control input σ k , k = 1, ..., 4 and Table 1 gives the optimal control parameters. The corresponding optimal objective function is 5.4791 × 10 −4 , which is very small and close to the optimum. In addition, we test the case of p = 1, which can be viewed as a constant value control strategy or an experimental experience strategy, and obtain the optimal control parameter 0.0249. The corresponding optimal objective function is 0.04852, which is much higher than the optimal objective value of the case of p = 4. Figure 3(a) shows the optimal colloid density distribution at t = 50s, t = 100s, t = 150s, t = 200s. Figure 3(b) shows the numerical errors between target colloid density distribution and optimal colloid density distribution along the space, which denote searching optimal control sequences approximate the optimum very well. Figures 4(a) and 4(b) show the spatiotemporal evolution of colloid transport by the optimal control input u(τ ) and the corresponding contour map, respectively. The colloids accumulate near the leading edge of the migrating colloidal front, which creats a pluglike colloidal "wave", and move toward the dead-end channel shown in Figures 4(a) and 4(b).
Note that the gradient-based optimization method of FMINCON in MATLAB and other nonlinear programming algorithms of advanced software are designed to find local optimum, and there is no way to guarantee global optimum by only applying such algorithms. Furthermore, the initial guess is usually related to the   Table 2 gives the optimal control parameters, which are close to the results of the case of p = 4.
6. Conclusion. In this paper, we have studied 1D colloid transport, which involves solute diffusion equation and the colloid transport model, and propose a PDE optimal control problem with the boundary control u(τ ) (solute concentration in main   Table 2. Optimal control parameters u(τ ) = σ k , k = 8 channel physically) during the colloid transport. A new effective computational method which involves control parameterization, variational method and numerical computation is proposed to design active optimal boundary control for this process. Then a gradient-based algorithm is proposed to solve the colloid transport PDE optimal control problem. The simulation results give the optimal control strategy and show the effectiveness and feasibility for the proposed control strategy. However, there are many challenging problems, such as model inaccuracy [31], computational physics [32], PDE-constrained optimization [9], topology design [15] in the microfluidics area. We plan to apply the iterative learning control scheme and linear quadratic (LQ) control scheme for the colloid transport PDE process [23]. For the time intervals [τ k−1 , τ k ), k = 1, ..., p, the time-scaling method which adaptively selects the optimal switching time instants can be applied [36,14]. Strictly speaking, since the colloid transport is a switching system with time delay, which involves the bubble case and solute input case, a computational optimal control problem of switching system with time delay could be considered [13,33]. Moreover, we are now building an experiment towards active control on microfluidics. The long term purpose of this paper is trying to implement the optimal control strategy for lab-on-chip colloid migration, and the experimental work is still in progress.
By using integration by parts, we can obtain Since C(x, 0) = Ψ(x), ∂C(1,τ ) ∂x = 0 and C(0, τ ) = u(τ ), (32) can be simplified to become To continue, we let µ : [0, 1] × [0, Z] → R denote another Lagrange multiplier function. Base on the dimensionless diffusion model (1) and the dimensionless colloid transport model (5), we have (34) Then, by using integration by parts again, we can obtain Since (35) can be simplified to give We now compute the gradient of J p at a fixed point σ = σ 1 , . . . , σ p ∈ R p , and let θ ∈ R p be an arbitrary vector as follows: The augmented objective function takes the following form: (38) Under the constraints of the solute diffusion equation (1) and colloid transport model (5), the optimal control problem to minimize the objective function (8) can be transferred into an unconstrained optimal control problem by introducing the Lagrange multiplier functions λ(x, τ ) and µ(x, τ ).
Thus, the perturbed augmented objective function (38) takes the following form: dxdτ.