When simulating large channel protein systems, lack of good initial values often causes convergence problems in the solution of the nonlinear steady-state Poisson-Nernst-Planck (PNP) model, a large scale coupled problem. Standard continuation techniques utilizing the time-dependent PNP solutions can address this convergence issue but they are often inefficient. We propose a reduced model describing a near- or partial-equilibrium state as an approximation of the original PNP system (describing a non-equilibrium process in general). Based on the reduced model, we design three initialization methods for the solution of the PNP equations under general conditions. These methods provide the initial guess of the PNP system by solving a specifically designed Poisson-Boltzmann- or Smoluchowski-Poisson-Boltzmann-like equation(s). Numerical experiments on 3D molecular systems (from protein databank) show that these methods can effectively reduce the number of Gummel iteration steps and/or the total CPU time when solving the PNP equations, and especially do not need continuation approaches anymore. Our numerical tests also show that as one of the initialization methods, the "linear approximation method" can even produce very close results such as current-voltage curves to that from the original PNP model when the membrane potential is not large.