## 1 Introduction  
In nature and engineering applications, there exist some quantities that potentially remain to be non-negative, such as the populations of organisms within ecological systems [1], the propagation rate of a signal in network communication [2], and the concentration of liquid in chemistry [3]. Such a system characterized by non-negative quantities is commonly known to be positive for which the state variables and output trajectories reside in the first quadrant whenever the inputs and initial states are non-negative [4]. The studies concerning positive systems are derived from [5], based on which significant results such as reachability and controllability [6], Kalman-Yakubovich-Popov [7] have been developed. Due to the positivity, the state variables remain in the non-negative cone instead of the entire linear space [8]. This property leads to the obstacle of applying the well-established general theories to positive systems, which promotes specialized research including state-feedback control [9], output-feedback control [10, 11], observer design [12, 13], etc. On the other hand, actual systems bear complex dynamics and the order of exact models is relatively high, which brings about great difficulties in the control and implementation of related systems. Naturally, the model reduction problem aiming at approximating a higher-order model by a lower-order one with a sufficiently small error has received considerable attention [14, 15, 16, 17, 18], and various model-reduction techniques have been established. The traditional balanced-truncation method was extended for the inhomogeneous initial condition case [19]. A lower-order stable transfer function was explored to approximate a given transfer function with the \(H_{\infty}\) norm [20]. A complete characterization of all optimal Hankel-norm approximations to a rational transfer function was provided in [21]. Nevertheless, two challenges arise when attempting to adapt the developed model reduction approaches to positive systems. On one hand, when approximating a specified positive system, we expect that the reduced-order model is preserved to be positive, which imposes additional constraints on system matrices. On the other hand, the approximation error is usually characterized by specific criteria of an associated error system, which can be transformed into a type of bilinear matrix inequalities (BMIs) involving the coupled terms between the Lyapunov variable and unknown system matrices [22]. An intuitive idea to resolve this issue is to transform the bilinear constraints into linear ones. In addition, the congruent transformation method [23] and the D-K iteration method [24] have been proposed to find feasible solutions in some cases. Nevertheless, the above results did not take the intrinsic positivity constraints into account, and solving the non-convex stability condition and the positivity constraint simultaneously remains a major challenge. This paper explores the positivity-preserving model reduction method for positive systems. A reduced-order positive model is sought to capture the dominant dynamics of the original high-order system within a prescribed error bound. The sufficient and necessary condition with the form of BMI for the positivity-preserving model reduction is developed in a way such that the approximation error is minimized within a guaranteed level. With the aid of the inner-approximation strategy, we develop a successive convex optimization (SCO) algorithm, under which the desired reduced-order model is achieved by sequentially solving the approximated convex problems. Overall, the main contributions of the paper are stated as follows:  
1. This paper proposes a positivity-preserving model reduction scheme, which guarantees the \(H_{\infty}\) performance of the resulted error system and preserves the positivity simultaneously. 2.