The production of polymeric materials occupies one of the leading positions in today’s chemical industry. The production of polymeric products relies on the process of repeated addition of monomer molecules to the active sites of a growing chain. If two or more monomers are used as starting compounds, the process is referred to as copolymerization. In particular, synthetic rubber is produced by this mechanism. In view of the fast development of computers and extensive industrial use of polymerization processes, the issues of their mathematical simulation are of current interest. Building a mathematical model allows one not only to predict the properties of a product but also optimize the production process [1].
In describing the mathematical model of the copolymerization processes of under the understanding of chemical characteristics of conventionally are two approaches: kinetic and statistical. Kinetic approach is classical in solving problems of chemical kinetics and has successfully established itself not only in the study of physical and chemical phenomena, but also in the optimization of technological processes in the chemical industry. The basis of this method is the preparation and the numerical solution of kinetic equations for the concentrations of all types of particles involved in the process. These equations are derived from the conditions of the material balance for each component reactions involving the law of mass action, which determines the rate of formation and disappearance of this component.
As the number of monomer molecules can reach tens of thousands, the kinetic scheme, which includes all the main reaction occurring in the system is reduced to an almost infinite system of differential equations. Directly to solve such a system is not possible, so when you write the model equations, it is converted to moments of the molecular weight distribution [2]. Using such a simplification of the system allows us to calculate the average molecular characteristics of the resulting product [4], as well as under conditions of uncertainty of the rate constants of some elementary reactions successfully formulate and solve the inverse problem [3]. However, to get the picture changes the molecular weight distribution and study the composition of the product obtained in this approach is no longer possible.
A statistical approach based on the fact that the set of polymer molecules in a given sample it explicitly assigned a statistically equivalent set of realizations of a random process. That is, each molecule explicitly or implicitly treated as a separate realization of specific random process conditioned movement along the polymer molecules, and the probability of this realization is considered appropriate to her equal share of all the other molecules in the reaction system.
For realization a simulation modeling approach copolymerization process we apply the method proposed in 1977 by the American physicist Gillespie [4]. Describe the algorithm of the following sequence of steps.
1. For build a model of butadiene-styrene copolymerization, let us assume that the reactivity of the active center at the end of a growing chain is determined by the nature of the terminal unit. Then the kinetic scheme of butadiene-styrene copolymerization can be described by the following steps:
initiator decay
initiation of active centers
chain growth
chain transfer
chain termination by disproportionation
chain termination by recombination
where 
; M1, M2 are the monomers of the first and second type; Pn,m and Qn,m are the active and inactive polymer chains with length m + n, comprising m units of the M1 monomer and n units of the M2 monomer, respectively; ki, kp, kreg, kd, kr are the reaction rate constants of initiation, growth, chain propagation, disproportionation, and recombination elementary stages, respectively; A(β) = {1, if β = 1; else 0}; B(β) = {1, if β = 2; else 0}.
2. Transform the experimental rate constant of elementary reactions to stochastic rate constants according to the following equations:
for first order reactions;
for bimolecular reactions between different species (V is the reaction volume, NA is the Avogadro’s number).
3. Then calculate the reaction rate for every reaction according to the equation:
(1)
where 
is the rate constant of the i-th reaction in which reagents A and B participate; XA, XB are the numbers of reagent’s molecules.
The total reaction rate is then calculated as the summation of the individual reaction rates:
(2)
where n is the number of elementary reactions forming kinetic scheme of the copolymerization process.
4. Then the probability of any reaction taking place at a given time is calculated by the following equation:
(3)
It is apparent that 
.
5. Generate a random number r uniformly distributed between 0 and 1 and pick up such value k that the inequality took place:
(4)
Consequently, reaction under an index k has to result from an imitating choice.
6. Continuing reasoning similarly we will build all scheme of carrying out reaction.
Algorithm has been tested for the process conditions emulsion copolymerization of butadiene and styrene. The obtained number-average and mass-average molecular masses showed a satisfactory agreement with the results from the kinetic model of the process.