Abstract

A wave-piercing catamaran (WPC) is a high-performance ship that has been developed in recent years. Compared to other common ships, the WPC has higher lateral stability, larger deck area, lower oil consumption, and higher speed. However, under rough seas and at high speeds, the coupled heave/pitch motions of the WPC can easily produce coupled oscillations that seriously affect its seaworthiness. To solve these problems, a ride control system was designed for the WPC in this study. This system comprises two T-foils, two trim tabs, and a catamaran motion controller. The H2/H∞ controller was designed based on memory-based particle swarm optimization to alleviate the coupled oscillation resulting from heave/pitch motions. Numerical simulations were conducted to validate the proposed method, and the results showed that the proposed motion controller could obviously improve the sea-keeping performance of a WPC.

1. Introduction

A wave-piercing catamaran (WPC) is a fast, high-performance ship that has been developed in recent years [1]. Its main advantage is that, in rough seas, it can pierce a wave rather than ride on it, thereby significantly reducing the wave-added resistance and speed loss. Therefore, a WPC normally has a maximum speed exceeding 40 knots [2]. Further, because a catamaran is much wider than a single-hull vessel, it has relatively much higher lateral stability. Finally, under adverse sea conditions, it shows ∼15% lower roll angle and can withstand large waves; these are useful characteristics in practice [3].

However, a WPC, owing to its higher speed, is more susceptible to disturbances caused by wind, wave, current, and other relative elements. In particular, the influence of heave and heave motions should be considered. Coupled heave/pitch motions have direct and negative effects on the seaworthiness of a ship. For example, harmonic shake can threaten the crew, cargo, and equipment on a ship. Therefore, studies have investigated the longitudinal fin performance and developed vertical stabilizers for high-speed catamarans. Further, T-foil and flaps have been used as subsidiary mechanisms in WPCs to improve their dynamic performance [4].

Active fin stabilizers are commonly used for ship roll motion control, and their control methods, including LQ control algorithm [5], model predictive control (MPC) [6], neural network and fuzzy logic algorithms [7], and robust control algorithm [8], have been studied extensively.

However, few studies have focused on coupled heave/pitch control. Since 2001, the Giron-Sierra team at the University of Madrid has been studying the hydrodynamic characteristics of the T-foil and trim tab for high-speed passenger ferries ferry through extensive pool experiments [9, 10]. A classic proportional-integral-derivative (PID) controller was used for the actuators [11, 12]; however, the joint control algorithm did not fully address coupling problems. Although the control effect was stable, the control system performance should be improved further using modern artificial intelligence control theory and corresponding strategies. A ride control system (RCS) is considered feasible for improving vertical motion. Further, studies have already designed an RCS [13] and a longitudinal motion controller [14] for a wave-piercing boat; therefore, an RCS can also be considered feasible for improving longitudinal motion. We also design a ride control system for WPC by NSGA-II [15] and use μ synthesis theory and improved GA for its controller design [16].

However, with rapid developments in ocean engineering, military science, and modern industry, ship RCSs are becoming more complex and system performance requirements are increasing. An important factor hindering the acquisition of high-performance systems is the uncertainty of systems and external disturbances. Uncertainty can cause deterioration of control quality and even system instability. Therefore, uncertainties and their impact on accuracy must be considered when designing control systems. In this light, robust control theory is used to solve ride control problems.

In robust control, the H2 and H∞ norms are the most widely studied measures of robust performance indicators. Hybrid H2/H∞ control considering a system’s H2 and H∞ performance indicators is a type of multiobjective control problem. However, owing to the nonconvexity of the conventional hybrid norm control problem, traditional mathematical methods cannot effectively solve this problem. Further, a multiobjective optimization problem does not have only one solution; instead, its solution is the set of all noninferior solutions or Pareto optimal solution.

In recent years, the linear matrix inequality (LMI) method has been used most commonly to solve the multiobjective H2/H∞ control problem. The LMI method requires the set of parameters satisfying constrained LMIs to be convex. Therefore, some additional constraints are inevitably introduced to solve the control problem, making the optimization results conservative. If the H2/H∞ hybrid multiobjective problem satisfies the time-domain dynamic performance index and has stable poles, the system can simultaneously have satisfactory robustness, steady-state performance, and dynamic performance. Therefore, the study of this problem has important application value; nonetheless, it remains a very difficult multiobjective control problem. Because this type of multiobjective control problem has both time-domain characteristics and frequency-domain characteristics of the optimal design of the multiobjective robust control system, it is difficult to solve using the LMI method. However, the robust control system design using a memory-based particle swarm optimization algorithm (MBPSO) can solve a multiobjective optimal control problem with both time- and frequency-domain characteristics, and because its constrained parameter set does not need to be a convex set, the conservativeness of the system is reduced. Therefore, we apply the multiobjective particle swarm optimization (MOPSO) algorithm to the design of a robust control system.

In this study, we establish an uncertain mathematical model for a WPC and then propose an MBPSO-based H2/H∞ controller for the RCS with a suitable T-type hydrofoil (T-foil) and trim tab. Finally, we conduct a numerical simulation to validate the proposed controller. The study results should have academic and industrial implications for improving the RCS of current catamarans.

2. Model of Coupled Heave/Pitch Motions of a Catamaran

Assuming a homogenous fluid medium in quiet seas, for small disturbances in the vertical plane, the catamaran’s frequency domain formulation can be expressed as [17]where , respectively, represent the displacement, velocity, and acceleration of the vessel; the indices k and j indicate the j-mode oscillatory motion caused by the k-direction force; M is the rigid body mass matrix; A and B are the added mass and damping matrix, and their values are based on the shape of the ship, forward speed, and water depth; C is the restoring matrix; and F is the excitation force.

For heave and pitch, the coupled two-degree-of-freedom (2-DOF) equations arewhere the numerals 3 and 5 denote the third and fifth DOFs that, respectively, indicate the heave and pitch movements. Correspondingly, F3 and F5, respectively, denote the heave force and pitch force of the target WPC. The hydrodynamic coefficients in these equations are mainly related to three factors: the static buoyancy effect of the ship, potential flow motion of the ship with a free water surface caused by a disturbance, and the effect of viscous flow on the naked hull.

The hydrodynamic coefficients related to the static buoyancy of ships can be obtained easily, whereas those caused by the potential flow are mainly calculated using the 2.5D potential flow theory [18, 19]. The calculations of the hydrodynamic coefficients related to the viscosity and the actuators are based on experimental and semiempirical data. Finally, the total hydrodynamic coefficient is obtained by adding each individual hydrodynamic coefficient described above. In the above equations, only C is a constant; A and B depend on the encounter frequency of the waves.

Figure 1 shows the catamaran model investigated in this study. Further, Table 1 lists the main characteristics of the ship. The hydrodynamic coefficients can be obtained using the 2.5D strip method. Figure 2 and Table 2, respectively, show the added mass and the damping and restoring coefficients of the WPC when its speed is 40 knots.

3. Model of T-Foil and Trim Tab

Wave-induced vertical accelerations may cause seasickness and serve as an important limit on the operability of high-speed vessels. However, motion control may be effective in reducing the heave and pitch of a high-speed vessel. T-foils and trim tabs are used as parts of control systems for multihull semidisplacement vessels. Trim tabs are located at the transom stern and T-foils, in the forward part of the vessel. These two parts are the two actuators of the RCS. Because vertical ship motions are large in the bow, it is advantageous for a heave and pitch damping device to be placed close to the bow. The structures of the actuators, namely, the T-foil located below the keel and the trim tab located behind the ship, are, respectively, shown in Figures 3 and 4.

The dynamic lifts and moments acting on the foils and flaps can be expressed aswhere LT and LF are, respectively, the mean coordinates of the T-foil and flap; CL is the lift coefficient; αT, αTF, and αF are, respectively, the attack angles of the T-foil, flap, and trim tab; AT is the area of each fin on the T-foil; AF is the area of each flap; ρ is the water density; and U is the speed. For a T-foil with a flap, when the foil shape is chosen, CL only relates to αT and αTF.

The principal parameters of the two actuators installed on the ship, namely, the T-foil and the trim tab, are, respectively, shown in Tables 3 and 4.

Figure 5 shows the corresponding relation between lift coefficients and attack angles in the towing experiment and CFD method.

4. Uncertain Model of Coupled Heave/Pitch Motion of a Catamaran with T-Foil and Flap

When a ship moves in the ocean, it is affected by uncertainties such as wave disturbances. Further, the dynamic performance is an important aspect to determine the vertical motion performance. Therefore, it is important to fully consider the robustness and dynamic performance of the system when synthesizing the control system. However, when a deterministic robust control method is used to synthesize uncertainties with certain dynamic performance requirements, the design process becomes complex and conservative and may sometimes have no solution. Therefore, we use the PSO method for robust control problems to control the vertical motion of ships. The close-loop control system is described aswhere Fh and Mp are, respectively, the control force and moment for heave/pitch motions. Figure 6 shows the T-foil and trim tab and the corresponding forces and moments to clarify the system and the relations among the parameters in these equations.

For heave/pitch motions, the space representation iswherewhere x is the state vector, y is the output vector, u is the input vector (total control forces and moments generated by T-foil and flap that need to be allocated and optimized), F is the disturbance vector, and E and A are coefficient matrices for the nominal system.

The vertical motions of a WPC in rough seas are affected by both wave disturbances and other uncertain factors such as the uncertainty of system parameters. Therefore, the system’s dynamic performance, especially its robustness performance, plays an important role.

In equation (7), A and B are mainly related to the potential flow and viscous flow of water, respectively; these are affected by many factors and are difficult to obtain accurately. C is mainly related to the ship structure, load variation, and water viscosity; these parameters also change constantly. The other parameters are related to the T-foil and trim tab. This part of the perturbation is mainly caused by the uncertainty of CL and speed U. Thus, the uncertainty model of the parameters can be obtained as follows:where

We can translate the real system to the following standard form of H robust control:where

The disturbance is , where

Then,

Let . Finally, the standard robust control form is obtained as follows:

We rearrange the equations to highlight the uncertain parts as follows:

This gives

The standard form of the H2/H∞ control system for the vertical motion of a catamaran is

The H performance of the vertical motion control system of the catamaran is mainly related to the pitching angle and heave displacement. Control effects and the robustness of the system to disturbances and perturbations must be improved. For the H2 performance of the system, the heave value, heave velocity, pitching angle, pitch velocity, and constraints on the control force/moment need to be considered. Accordingly, the performance evaluation signal for the ship’s longitudinal motion control system is defined as

The problem in H2/H control is to design a controller K to make the closed-loop system asymptotically stable. Further, the H norm of the closed-loop transfer function from ω to z should not exceed a given upper bound γ1 to simultaneously ensure the robust stability of the closed-loop system to uncertainties and minimize the H2 norm of the closed-loop transfer function from ω to z2:

Then, the state feedback regulator K has the following structure:

The state feedback regulator coefficients for the T-foil are k11, k12, k13, and k14 and those for the trim tab are k21, k22, k23, and k24. These coefficients depend on the heave/pitch motions, heave/pitch velocities, and heave/pitch acceleration. Increasing the control coefficients k11 and k21 reduces the heave motion period of the ship. Increasing the control coefficients k12 and k22 reduces the heave response amplitude operator (RAO) of the ship. Increasing the control coefficients k13 and k23 reduces the pitch motion period of the ship. Finally, increasing the control coefficients k14 and k24 reduces the pitch RAO of the ship, that is, it reduces the ship’s vertical metacentric height.

5. State Feedback Regulator K Optimized by MBPSO

Traditional methods like LMI are usually known to be conservative in the design of multiobjective robust control systems. Therefore, PSO is applied to the design of the multiobjective robust control system. Toward this end, this study proposes a PSO-based H2/H∞ robust control design method. The simulation results show that compared with the traditional LMI method, the robust controller obtained using the MOPSO algorithm has better robustness performance index and reduces the conservativeness of the robust control design.

5.1. Objective Function

The objective of PSO is to update the state feedback regulator to obtain control signals (αT, αF) from the H2/H∞ controller. By assigning corresponding weights to the system’s H2 and H∞ performance, the multiobjective optimization control problem is transformed into a single-objective optimization control problem that needs to be solved. For the antivertical effect of a WPC, this study considers the minimum vertical motion as the optimization objective. The corresponding objective function can be expressed as follows:

Further, the constraint conditions are

The first and second terms in the objective function, respectively, represent the H∞ and H2 norms. Parameter optimization can be performed using various methods such as the simplex method and the gradient method. However, the former can easily fall into a local optimal solution and the latter can only be used with a continuous-differential performance index function, making it unsuitable for our case. Compared with other intelligent optimization algorithms, the PSO algorithm affords many unique advantages. First, the principle and structure of the PSO algorithm are much simpler and easier to implement than other algorithms. Second, the PSO algorithm model contains relatively few parameters that can be adjusted conveniently. Third, one of the most important characteristics of the PSO algorithm is that it can memorize and share location information for the population, and therefore, an individual and a global optimum can be obtained according to the position and speed. Finally, the PSO algorithm enables multiple particles to simultaneously search for the optimal solution in a complex solution space, thereby affording better parallelism.

5.2. MBPSO

PSO is an evolutionary computation method to simulate social behavior that has attracted increasing attention in recent years [20]. PSO performs better than a genetic algorithm (GA) in solving continuous optimization problems. The main problem PSO faces in a dynamic environment, like other evolutionary computation methods such as GA, is that all particle swarms tend to converge gradually with the iteration of the algorithm, making it impossible to respond quickly to changes in the environment.

Many studies have shown that when the environment does not change very drastically, transmitting evolutionary information from the previous population to the new population can effectively improve the performance of the algorithm and enable adapting to the changing environment more quickly. Memory strategies, especially explicit memory methods, have been widely used in GAs for solving dynamic optimization problems. These methods can guide the evolution of new populations by storing the best solutions obtained in the running process and reusing them when necessary. In this study, an explicit memory strategy is introduced into PSO to improve its ability to solve dynamic optimization problems. In fact, PSO itself has a certain memory function. According to formulas (23), each particle iteratively updates the best historical solutions found by itself and by all “partners” in its neighborhood during PSO operation. The PSO with the introduced explicit memory strategy is called the memory-enhanced particle swarm optimization (MEPSO) algorithm:where represents the number of particles in space; represents the dimension of the solution space; acceleration factors are positive constants that are set as ; and are random numbers between 0 and 1. Further, when , then , and when , then , where , respectively, represent the minimum and maximum velocities of particles in flight. In addition, when , then , and when , then , where , respectively, represent the nearest and farthest position ranges of particles in space. Finally, is the inertia weight, and its magnitude reflects the influence of the particle’s historical velocity information on the current velocity.

Intuitively, when an optimal solution reappears at or near the locations it has visited, the memory can help recall these locations and make the population move quickly to the new optimal solution. The memory can also help maintain the diversity of populations and guide the rapid evolution of reinitialized populations into viable regions. However, the introduction of the memory strategy will also cause the algorithm to attach too much importance to the development of knowledge acquired in the past while ignoring its search ability; this may mislead the evolution of the population and hinder its exploration of new search spaces and discovery of new peaks. Allocation optimization is a dynamic optimization problem, and standard PSO may not find the best solution. Therefore, we combine the advantages of the complete initialization method and the memory method and propose a memory-enhanced three-island PSO algorithm for solving dynamic optimization problems [21, 22].

The island model was originally designed for parallel evolutionary computation. Each subpopulation was isolated from each other and evolved concurrently, and the best individuals were migrated to adjacent subpopulations according to a certain period. The proposed PSO algorithm is a “three-island model” in which the whole particle swarm is divided into three parts: memory particle swarm, searching particle swarm, and developing particle swarm shown in Figure 7. The memory particle swarm is used to record the best solution that all particles can find in a certain period of time. When the environment changes, especially when extreme points return to these positions, this swarm can help the population quickly find these new extreme points. This swarm also monitors environmental changes. If a particle’s fitness changes and then when the external environment changes, all particles need to be revalued. The search particle swarm is used to ensure the population’s ability to explore space; it passes its best particles to the memory particle swarm every time for a period of time and then reinitializes. The developing particle swarm obtains information from the memory particle swarm, introduces some particles to the memory particle swarm when it is initialized, and periodically transfers its best particles to the memory particle swarm; further, when the environment changes, it needs to be reinitialized.

MBPSO is used for the H2/H∞ controller for the T-foil and trim tab. Generally, the solution is to find the inputs of multiple actuators on the premise that the constraints and the optimization objectives are met simultaneously. The PSO is used to update the state feedback regulator that depends on the objective function given in equation (20) with constraints. Therefore, for the studied H2/H∞ controller problem for the T-foil and trim tab, the position information of each particle is the size of (k11, k12, …, k24). In other words, the velocity information of each particle is the change rate of (k11, k12, … , k24). The function value of each particle is the objective function of the control system. The optimal allocation of coordinated control using the MBPSO algorithm involves the following steps:(1)Initialize algorithm parameters, including population size, solution dimension, and maximum number of iterations(2)Initialize the position and velocity of all particles for the problem of coordinated control allocation of the RCS, that is, initialize the positions and velocities of K’s each parameter(3)If the fitness of one particle changes, then develop a particle swarm, search for all particles in the particle swarm that need to be reevaluated and regenerate the exploited particle swarm(4)Update the search particle swarm and developing particle swarm(5)If the memory is updated, the best particles obtained by the searching and developing particle swarms are passed to the memory particle swarm and go to step 6; otherwise, go to step 7(6)Update the memory particle swarm and reinitialize the search particle swarm according to equation (23)(7)If the maximum number of iterations is not reached, repeat the calculation in step 2; otherwise, stop the calculation, as the global optimal solution of the objective function is obtained.

6. Simulation and Results

In this study, MATLAB Simulink was used to simulate the process and validate the proposed control strategy [19]. The model is divided into four parts: WPC vertical motion model, robust controller model, MBPSO model, and T-foil/trim tab model. Figure 8 shows the control block diagram.

6.1. Determination of Simulation Conditions and Related Parameters

For the catamaran with a T-foil and trim tab as described before, the simulation conditions and related parameters are as follows: ship speed, 40 knots; heading angle, 180; and sea state, 5. When the speed of ship is 40 knots, the hydrodynamic coefficients are different under different encounter frequencies. Considering the simplicity of calculation, the amount of calculation is small, and the accuracy can also meet the needs of modeling. In this paper, the method of dominant frequency point is mainly used in the modeling; considering the influence of encounter frequency on frequency response, the dominant frequency point is chosen as 1.26, because the frequency responses of pitch and heave are the biggest, and the model is as follows:

Further, waves are generated using the PM spectrum; Table 5 shows the wave parameters.

We use PD control and robust control with LMI, PSO, and MBPSO in the simulations. The initial position of each particle in the MBPSO algorithm, that is, the initial attack angles of the T-foil and trim tab, is . The maximum number of iterations of the PSO is set to , population size is set to N = 40, and c1 = c2 = 2. The solutions’ dimension is D = 8 for each coefficient. We set the same wave-induced forces/moments and system parameter perturbations in all three methods. In MBPSO, the number of update iterations plays a very important role in optimization. Several simulations indicate that the vertical motion reduction is higher with more frequent updates. However, the vertical acceleration exceeds the maximum value for an overly large number of update iterations. Therefore, the number of update iterations was set in the range of 1–25 in the simulations. The exploit-swarm and explore-swarm had 20 particles, and the memory particle swarm had 10 particles.

6.2. Model of Irregular Wave Disturbance Force/Moment and Performance Evaluation

Generally, irregular waves can be considered to comprise an infinite number of constituent waves having different amplitudes, frequencies, and directions and the characteristics of phase disorder. These waves are combined to form the wave spectrum that describes the distribution of wave energy relative to each component wave; therefore, it is also called the “energy spectrum.” In this study, the Pierson–Moscowitz (P-M) spectrum is used to simulate wave motion:where the unit of S (ω) is m2s and is a significant wave height.

The force/moment acting on a ship with different wave frequencies can be obtained from the added mass and the damping and restoring coefficients in Table 2 and Figure 2. We decompose the forces and moments of waves acting on the hull into the superposition of the disturbing forces generated by regular sinusoidal waves as follows:

The disturbance force at each frequency can be obtained from Figure 2. After adding the disturbance force/moment at each wave frequency, the force and moment of the wave acted on the ship are shown in Figure 9.

In this paper, two indexes are added to evaluate the control effect: the total percentage damping ratio, R and the statistical vertical motion reduction constant, S. R describes the overall motion decrease between controlled motion and uncontrolled motion:where subscript p represents pitch motion and subscript h represents heave motion.

In this calculation, the total vertical motion reduction ratios are evaluated as total pitch angles and heave displacements; however, it only focuses on total pitch angles/heave displacements reduction and does not make any distinction between small and large vertical motions in terms of proportions. The statistical vertical motion reduction constant (S) allows one to distinguish between small and large vertical motions. The higher the constant F is, the better performance the ride controller has. F can be calculated as follows:where mocp is the mean square value of controlled pitch angles, moch is the mean square value of controlled heave displacement, moup is the mean square value of uncontrolled pitch angles, and mouh is the mean square value of uncontrolled heave displacement.

6.3. Parameter Optimization Performance of MBPSO

Figure 10 shows comparisons of the heave and pitch between different parameter optimization methods. The maximum heave/pitch motion for PD control is 0.6 m/2.5 degrees, which is much larger than those of intelligent algorithms. PSO and MEPSO have the best optimization performance, and their heave/pitch motions are both within 0.3 m/1 degree. The maximum heave/pitch motion for LMI is 0.4 m/2 degrees; this is slightly worse than those of PSO and MEPSO.

From Table 6, we can see that MBPSO-based H2/H shows the best performance with biggest R and S and PD control shows the worst performance with lowest R and S.

Figure 11 shows a comparison of the fitness values between LMI, PSO, and MEPSO. This figure shows that MEPSO has the fastest convergence rate and the highest convergence precision. The improved optimization performance of MEPSO, compared to that of PSO, is attributable to the introduction of the three-island model that accelerates the local search and avoids falling into a local optimum solution. The MBPSO algorithm is updated when a new environmental condition appears, or the algorithm randomly chooses an iteration step to update. Therefore, some particles are randomly initialized to explore new solutions; further, the remaining particles are recalled from the memory to exploit past solutions.

6.4. Performance of H2/H∞ Antivertical Motion System

To verify the robustness of the proposed MBPSO-based H2/H controller, Figure 12 and Table 7 show the WPC’s heave/pitch motion with both the MBPSO-based robust controller and the conventional robust controller, where (i = 1, … , 12) is varied by ±30%.

By contrast, the heave/pitch motion of the conventional PD controller is within 0.8 m/3.5 degrees when is varied by 30%; this is much larger than that of the MBPSO-based H2/H controller. The two PSO-based H2/H controllers’ output signals almost overlap with each other, and both of them have nearly minimum heave/pitch motion within 0.5 m/1 degree. Further, the heave/pitch motion of the LMI-based H2/H controller is within 0.8 m/2 degrees; this is much larger than that of the MBPSO-/PSO-based H2/H controllers. MBPSO also has the largest R and S. This indicates that the MBPSO/PSO-based H2/H controllers have much better robustness than the LMI-based H2/H controller.

6.5. Other Analysis Results

To prove that the proposed controller is feasible under other sea states, we simulate two other sea states. Table 8 lists the parameters for these sea states. For simplification, only PD control and the MEPSO-based H2/H controller were simulated. Figure 13 and Table 9 show the results.

According to the simulation results presented above, when the MEPSO-based H2/H controller and PD controller were added, the wave-induced heave/pitch motion decreased by 90%/65% on average with MEPSO-based H2/H controller, indicating that this controller was theoretically feasible and was much better than PD controller.

7. Conclusion

In this study, the WPC’s uncertain vertical motion model is established and a calculation method for the lift force and moment produced by the T-foil and trim tab is provided. These are substituted into the WPC’s longitudinal motion differential equations for the design purpose of making the restoring force and torque generated by the T-foil and flap offset, and the wave force and torque are generated by waves. Accordingly, the uncertain longitudinal motion equation of a WPC with an RCS is established.

After the motion equation is established, the controller is designed. This study discusses the origin, development, and basic details of H2/H control theory and PSO. When designing H2/H controllers, the most important factor is the selection of the H2/H state feedback regulator. Then, parameter optimization based on the controller for the T-foil and trim tab is performed using MBPSO. The RCS reduced most of the wave-induced force and torque. Compared to PD control and LMI-based H2/H control, the proposed MBPSO/PSO-based controllers provide much better antivertical motion effects. This study provides an important reference for the control system of a WPC. In the future, we will explore more effective methods to facilitate the controller and use simulations to verify the design. Further, the above calculations can also be applied to the design of other types of controllers.

Data Availability

The data in this work were obtained by means of MATLAB Simulink tool, and the same data are believed to be reproduced on employing the same parameters and tool and similar procedures mentioned in the paper.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This project was supported by the Natural Science Foundation of China (grant no. 61603110).