With the growing demand for efficient co-utilization of renewable and traditional energy sources, co-pyrolysis of biomass and coal has gained significant attention as a green energy conversion pathway. However, the process involves complex multivariate interactions, and traditional modeling methods struggle to identify key variables and their interactions accurately. This study aims to develop a multivariate modeling and collaborative optimization framework for co-pyrolysis to enhance product distribution prediction and process control efficiency. Based on experimental data from various biomass and low-to-medium rank coal co-pyrolysis, a multiple linear regression model incorporating main effects and interaction terms is constructed, with hexane insolubles introduced as a structural product indicator. Significant variables and synergistic effects are analyzed using methods such as Kendall Tau rank correlation and Mann-Whitney U tests. Further, Pareto frontier analysis and Monte Carlo simulations achieve dual-objective optimization of tar yield and hexane insolubles. The results indicate that there is a significant interaction between hexane insolubles and the blending ratio, enabling optimal coordination of tar yield and product structure. Specifically, hexane insolubles show a strong positive correlation with char yield, while their correlation with water and char is relatively weak. Pareto frontier analysis revealed that at an optimal biomass blending ratio of 20%–30%, the tar yield increased by approximately 12% compared with coal-only pyrolysis. Furthermore, Monte Carlo validation based on 100 randomly generated feasible solutions demonstrated that the Paretooptimal set outperformed random solutions in all 100/100 trials, confirming the robustness and generalization ability of the model. These findings demonstrate that the proposed framework not only achieves high predictive accuracy but also provides practical guidance for optimizing co-pyrolysis systems.