Научная статья на тему 'Optimization of the shape of the Pareto set in the problems of multi-criterial programming'

Optimization of the shape of the Pareto set in the problems of multi-criterial programming Текст научной статьи по специальности «Математика»

CC BY
94
31
i Надоели баннеры? Вы всегда можете отключить рекламу.
Область наук
Ключевые слова
MULTI-CRITERIAL PARAMETRICAL PROGRAMMING TASKS / SET OF PARETO / METHOD OF SMOOTH PENALTY FUNCTIONS / OPTIMIZATION PROBLEM IN TERMS OF PARAMETERS / ЗАДАЧА МНОГОКРИТЕРИАЛЬНОГО ПАРАМЕТРИЧЕСКОГО ПРОГРАММИРОВАНИЯ / МНОЖЕСТВО ПАРЕТО / МЕТОД ГЛАДКИХ ШТРАФНЫХ ФУНКЦИЙ / ЗАДАЧА ОПТИМИЗАЦИИ ПО ПАРАМЕТРАМ

Аннотация научной статьи по математике, автор научной работы — Katendi Bula Axell, Umnov E. A., Umnov A. E.

В работе рассматривается схема использования метода гладких штрафных функций для исследования зависимости решений задач многокритериальной оптимизации от параметров. Приводится описание алгоритмов, основанных на методе гладких штрафных функций, решения задачи оптимизации по параметрам уровня согласованности целевых функций и выбора соответствующей формы множества Парето.I

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

n this paper, a scheme for using the method of smooth penalty functions for the dependence of solutions of multi-criterial optimization problems on parameters is being considered. In particular, algorithms based on the method of smooth penalty functions are given to solve problems of optimization by the parameters of the level of consistency of the objective functions and to find the corresponding shape of the Pareto’s set.

Текст научной работы на тему «Optimization of the shape of the Pareto set in the problems of multi-criterial programming»

Review of Business and Economics Studies 2018, Vol. 6, No. 1, 5-16

DOi: 10.26794/2308-944X-2018-6-1-5-16

Optimization of the Shape of the Pareto Set in the Problems of Multi-criterial

Programming

Bula Katendi Axel1, E. A. Umnov2, A.E. Umnov2 Moscow institute of Physics and Technology (State University) [email protected], [email protected]

Abstract

in this paper, a scheme for using the method of smooth penalty functions for the dependence of solutions of multi-criterial optimization problems on parameters is being considered. in particular, algorithms based on the method of smooth penalty functions are given to solve problems of optimization by the parameters of the level of consistency of the objective functions and to find the corresponding shape of the Pareto's set. Keywords: multi-criterial parametrical programming tasks; set of Pareto; method of smooth penalty functions; optimization problem in terms of parameters. JEL classification: C61, C65

I

n mathematical modeling, it is often necessary to formalize preferences for states of the modeled object that generates several independent target functions. According to the historically established tradition, in this case, it is customary to talk about of multi-criterial optimization problems. A finite-dimensional multi-criterial model is a mathematical model with N objective functions:

fk (x,umax k = [1,N], (1)

subject to maximization possessing at interior points of the set of elements x e En, and satisfying the following conditions:

y (x,u )< 0 i = [1, m ] (2)

where u e© c Er — vector of parameters of the model. It is assumed that the functions fk (x,u) and yt (x,u) are sufficiently smooth, i.e. they have continuous derivatives of a desired order in all their arguments.

The incorrectness in the general case of such a statement is obvious, since the element x that is extremal for one of objective functions, in general, is not such for others.

However, useful information can be obtained by successively solving the following problems with a criterion for finding an extremum on the set (2) of each of the functions (1) separately for k = [1, N ]:

fk (x,u)

^ max

x

subject to y¡ (x,u)< 0 i = [1,m] (3)

The objective function fk (x,u) is called improvable in the feasible point x0 (i.e., satisfying the condition (3)) if there is another feasible point xx, for which fk (xvu) > fk (x0,u).

It is clear that the solution of problem (3) for any k = [1, N] is un-improvable, or "non-ideal" at the point of view of the other objective function fk (x,u).

The concept of improving the multi-criterial objective function allows the feasible points to be divided into two subsets: for the first, all feasible points improve all objective functions and for the second, there are points for which the improvement of one function causes the deterioration of at least one other function.

The second subset is called a Pareto-type set or, simply, a Pareto set.

A general universal approach to the solution of multi-criterial optimization problems has not been proposed yet, but numerous approaches have been developed (see Fiacco & McCormick, 1968; Lotov & Pospelov, 2008), which limit the number of solutions.

For example, in the practical use of multi-criterial mathematical models, the set of independent objective functions is often replaced by a single one, thus passing to the standard problem of mathematical programming, allowing finding consistent or compromising solutions on the Pareto set in a certain sense.

Statement of the problem

In this article, the problem of finding an element on the set (2) that minimizes the gap between the objective functions will be considered as a compromise. In other words, this is a mathematical programming problem of the following form:

minimize p, subject to p > 0,

y (x,u )< 0 i = [l,m ], (4)

fk (x,u)> F* (u)-p k = [1,N]

whose solution will be denoted by p** (u) and x** (u) . Here Fk (u) = fk (x* (u),u) and p is called "mismatch value".

The problem (4) is naturally called a two-level parametric problem since in its formulation it contains solutions F* (u) k = [1, N] of problems with single criterion (3), which we call first level problems. In this case, both in the problems of the first and the second level, it is assumed that the vector of the parameters u e© is fixed.

It is clear that the extreme value of the mismatch between the criteria in the general case is determined by the properties of the Pareto set and depends on the parameters vector u . Therefore, it is natural to indicate the third level optimization problem for the models (1)-(2) as follows:

optimize the expression p** (u) subject to u e© (5)

which solution will be the vector of parameters u*** e© and the number p*** = p** (u***). In the present paper, possible solutions to problem (5) will be considered.

Solution method

Let us consider the problem of finding in the parameter space a standard method (for example, gradient) of finding the extremum of the mismatch value of the objective functions of the multi-criterial model (3)-(4)-(5).

The specificity of this problem is based on the fact that the formulation of the problem (5) (the upper level or third level) includes the dependence p** (u), the solution of the problem (4) (the second level) which in turn, depends on F** (u) = fk (x* (u),u) Vk = [1,N] — the solutions of the problem (3) (the lower level or first level).

The functions p** (u) and Fk (u) Vk = [1,N] in the general case (even for smooth functions fk (x,u) and y (x,u)) may be not differentiable, that why the use of any numerical method based on Taylor approximations is not possible.

It is proposed the use of the method of smooth penalty function to overcome this difficulty (see Umnov, 1975) and obtain a sufficiently smooth approximation dependences of p** (u ) and F* (u) Vk = [1, N].

It's assumed that the penalty function P(t,s), which penalizes the restriction s < C , satisfies the following conditions:

1 Vt > 0 and Vs , the function P (t, s) has continuous derivatives with respect to all its arguments up to the second order

2 Vt > 0 and Vs ,

dP . d2P — >0 ; —->0 ds ds2

(6)

3 P(t,s)> 0Vs and Vt> 0, and,

/ \ [+œ, s > 0 lim P (t, s) = <! t^+0 v ; I 0 , s < 0

(7)

When solving the third-level problem by an iterative method, for each step of the method, it is necessary preliminary to solve the problems of the second and first levels for a fixed vector of parameters u. Let us first consider a possible scheme for solving first-level problems. In fact, we will use an auxiliary function for the one-criterion problems (3), as follows:

Ak (t,x,u) = fk (x,u)-^P(t,yi (x,u)) Vk e [1,N]

(8)

i=1

while a sufficiently smooth penalty function P (t, s) satisfies conditions (6) and (7).

As shown in Zhadan (2014), instead of the smooth approximations x* (u) solutions of each task of problem (3), we can take xk (u) stationary points of the auxiliary function (8), defined like:

dA*

= 0 Vj e [1, n ]

(9)

or

f _ » dp ôyL = 0 dxj tï fyi dxj

Vj e [1, n ]

Since the condition of the second-level problem (4) includes the dependencies K (u) = fk (xk (u),u) V* = [1,N] which are not differentiable functions for aH their arguments, then for these dependencies it is also necessary to choose a smoothed approximation.

As an approximation, the auxiliary function calculated at a stationary point Fk (u) = Ak (t,xk (u),u) can be used, because (due to the properties of the penalty function method) its value for small positive t is close to the optimal value of the objective function of the k -th problem (3).

Standard optimization methods used for lower-level tasks, based on the use of continuous gradients or other differential characteristics, suggest that in addition to the solving system (9), these characteristics themselves can be found.

Let us demonstrate this using the example of calculating the derivatives of the function Fk (u) with respect to the components of the vector u of parameters.

As Fk (u) = Ak (t,xk (u),u), then according to the rule for differentiating a composite functionof several variables, we have:

ÔFl =3Al " dXj_ dup dup dXj dup

Vp e[l, r ]

Using (9), we have:

dFk ôAk t _ , s. \

-r^ = 1rjL(tXk (u),u) dup dup

Vp e [l, r ]

(10)

Note that the last simplification would be impossible if for Fk (u) a more natural approximation fk (x* (u),u) is used instead of the smoothing approximation Ak (t,xk (u),u).

Let us now look into the solution to the second-level problem. To make application of the penalty function method more convenient, the problem (4) is expressed as: maximize -p, subject to -p > 0,

yt (x,u)< 0

i = [l,m],

(11)

Yk (p,x,u)< 0 k = [1,N]

where Yk (p,x,u) = Fk (u)-p-fk (x,u)

The solution to this problem will be denoted by p** (u) and x** (u). Let us define the auxiliary function for the problem (10) as follows:

N m

e (- ^x,u )=-p- p (т, -p)-^p (т, Yk (p, x,u ))-YP (т, y>- (x, u))

(12)

k=1

i=1

replacing previously in Yk (p, x, u) the dependency Fk (u) by its smoothed approximation Fk (u).

For the set of variables {p, x1, x2,...,xn} , the conditions for the stationarity of the auxiliary function (12) will be:

dE dP N dP „

— = -1 + — + > -= 0

dp dp ¿=1

dE * ep_ f dPdy±. = 0

dx} kk dYk dx} tt dyt dx}

(13)

Vj = [1, N ]

Let the solutions of system (13) be p(u) and x (u), then, as a smoothed approximation of the dependency p** (u), we can use the function E (u) = -E (p, p(u), x (u ),u). The derivatives of this function by the components of the vector u ofparameters and the rule for differentiating a composite function of several variables give us:

dE dE " dE 3x, dE dp w ri -,

-=-+>--- +--— Vp e [1,r ].

dup dup dx, dup dp dup

From (13) we know that — = 0 and — = o Vj = [1, N ].

dp dx,

Then the last expression can be written simply:

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

If = lf(p,P(u),x (u),u) VP e[l,r] (14)

p p

Finally, we obtain formulas for the gradient components of E (u) in terms of the functions usedin the formulation of the multicriterial model (4)-(5) and the method of smooth penalty functions. From (12) it is obtained:

dE =JL dp_ dYk_ ap dyj_ dup k=i dYk dup k=i dup

d y dF df d F

where —^ = —L _J±. and the value of_k can be found from (10).

dup dup dup dup

Formulas (14) allow us to solve the third-level problem by applying any of the first-order methods, for example, conjugate directions. Note that second-order methods should also be considered here. However, this will be done at the end of the article, while now let us illustrate an example.

Proposed method in use

Let us consider multi-criterial mathematical model in which x = || x1 x2x3 ||re E3 is a vector of independent variables and u = || u1u2 ||re E2 is a vector of parameters. The problem is to maximize for x and u e© the functions:

U

xx

Figure 1. Geometric interpretation of models (6)-(7).

fi (x,u) = xi, f, (x,u) = x2, f2 (x,u) = x3 subject to x1 > 0, x2 > 0, x3 > 0, a1 (u1;u2) x1 + a2 (u1;u2)x2 + a3 (u^u2) x3 < b (u1;u2)

and where the functions a1 (u1;u2 ),a2 (u1;u2 ),a3 (u1;u2) andb (u^u2) are given by the condition below.

A valid region of the model (with an allowable fixed u) is a rectangular pyramid OABC. The Pareto set coincides with the face of ABC or is a part of it. A = || u1001|7 and B = || 0u201|7.

We assume that the set © in the parameter space is given by the condition that the sum of the lengths of the segments OA, OB and OC is constant and equal 3.

Applying the standard methods of analytic geometry, we find that for the compatibility of the system of model constraints, the existence of r > 0 is necessary such that:

a1 (ux,u2 ) = u2r a2 (ux,u2) = u1r a3 (u,u2) = u1u2 b(u,u2) = uxu2r

By choosing r such that r = 3 - u1 - u2, we assure that the set © will not be empty when 0.1 < u1 < 2.5 and 0.1 < u2 < 2.5.

The minimum value of the discrepancy between the criteria in this example depends on the form of the Pareto set, which is the triangle ABC, or a part of it. A graphic representation of the dependence of the error value of the objective functions on the parameters u1 and u2 is shown in pictures 2 and 3. Let us see with more details the properties of this dependence.

Clearly, the solutions of the first-level problems (3) for fixed u1 and u2 are:

f (x* (u)) = ux, f (x* (u)) = «2, f (x* (u)) = 3 -ux -«2.

Consequently, the task of the second level (4) — minimizing the discrepancy of the criteria, will have the form:

minimize p according to {x1, x2, x3,p} subject to p > 0 ,

x1 > 0, x2 > 0, x3 > 0, u2rxx + uxrx2 + u1u2x3 < u1u2r ,

x1 >u -p, x2 > u2 - p, x3 > r - p '

r = 3 - u1 - u2

His solution will be designed by p** (u1,u2).

Finally, the task of the third level (5) for our example will be: minimize p** (u1,u2) by {u1,u2} when 0.1 < u1 < 2.5 and 0.1 < u2 < 2.5.

It is known from the theory of mathematical programming that the properties of the dependence p** (u1,u2) are primarily determined by how the set of constraints of a model of the «inequality» type is divided into active and inactive ones, that is, the first of which are satisfied as equalities, and the second — as strict inequalities.

This separation depends on the values of the parameters of the model and its optimal variant determines the solution of the second-level problem.

First, suppose that the values of the model parameters initiate a conflict of all three criteria simultaneously. In other words, the improvement of the value of any one of the objective functions of the model is possible only if the values of all the others deteriorate.

In this case, the last five constraints of the second-level problem must be active, and we obtain

Figure 2. The quantity |OA| + |OB| + |OC is constant (graphic - 3D).

the following system of equations, which allows us to find the analytical form of the dependence p** («i,«2).

U2 rXi + u r%2 + Ui U2 X3 — Ui U2 r

x1 — u1 - p < X2 — u2 -p X3 — r-p

r — 3 - u - u2

Solving this system above, we have the analytic form of p which depends on u1 and u2:

2

P** (ui,u2 ) — -J-j"

u1 u2 3 - u1 - u2

The stationary points of p** (u1,u2) are || 1 1 ||r, || -3 3 ,y 3 -3 ||r and y 3 3 can be easily found, and for the first point the function has a local maximum with the value 2/3 according to the Sylvester criterion, not for the others because they don't satisfy the condition of non-negativity of the variables x1, x2, x3 and r .

The formula obtained is valid only in a certain area contained in © . An analysis of the isoline system shown in Picture 3, allows selecting five areas with different sets of active restrictions. Light lines determine the boundaries between the areas. The formula obtained above is valid only in area 4. In this area, the Pareto set of the model under consideration consists of the interior points of the triangle ABC.

Outside area 4, the formula for p** (u1,u2) is different. For the area 1, for example, p** (u1,u2) is found from the system of equations:

u rxj + u*i rx^ + u>i u ^3 — u>i u r x1 — 0

x2 — u2 -p x3 — r-p

r — 3 - u1 - u2

1

jti v ■VI №

J-1 T

J>r1 -j

("Si jti J5

j*I a

j 1-1.a

22

Jti u ■H H jti m jt*s

Jti u jtt-vj

Jti.M

Jti 99 MS

jts.ji

Jti 36 ■M E

(Ml

Figure 3. The quantity |OA| + |O5| + |OC is constant (isoLine view).

since the set of active restrictions in it is different: it contains the condition x1 = 0 , instead of x1 =u1 - p . We can easily prove that in the area 1

** / \ 2

p (u1,u2 ) = 1-1-

—+-

u2 3 - uy - u2

There are not stationary points for this dependency.

For areas 2 and 3, the arguments and results are similar. The Pareto sets in areas 1, 2, and 3 are the sides of the triangle ABC: BC, AC, and AB, respectively. Finally, we note that in area 5 the system of conditions (2) is contradictory.

In this case, the exact solution to the problem of the upper (third) level has the form:

2

Uy** = 1 , u2* = 1 , p** = 3

Let us now describe the method for solving the third-level problem for the variant of the multi-criterial model. We have:

objective functions: maximize by x = || xyx2x3 ||re E3

fy (x,u) = Xy f2 (x,u) = X2 f3 (X,u ) = X3

subject to xy > 0, x2 > 0, x3 > 0, ay (u) xy + a2 (u) x2 + a3 (u) x3 < b (u) where

ay (u) = u2 (3 - u - u2) a2 (u) = u (3 - uy - u2) a3 (u) = uyu2 b (u) = uyu2 (3 - u - u2)

with

u = yulu2 II e0 = ^«

0.1 < u < 2.5 J 0.1 < u < 2.51

Let us introduce the notations:

y = - x ^2 = - X2 ^3 = - X1

y4 = a1x1 + a2 x2 + a3 x3 - b and as a penalty function, we take P (t, s) = t exp ^ T j.

Then, the auxiliary functions for the one-criterion problems (8) will be:

k = 1,2,3

4 (y \

Ak (t,x,u) = Xk -£texp

j=i

v t /

(15)

The stationarity conditions of the auxiliary functions by the components of x give the following equations:

dAk R

âX7 = 8* +exp

y

v t /

- ajex^p\y4 I Vk = 1,2,3 Vj = 1,2,3

Finally, the derivatives of the smoothed approximating functions Fk (u) over the components of the parameter vector u are found from the relations (10)

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

— dA

( Fk L =iâ (t, xk(u ),u )■

du

da1 _ da2 _ da3 db ^

xk1+ xk 2 T + xk 3 t ~

dup dup dup dup

v p p P P /

exp| — ¡Vk = 1,2,3 Vj = 1,2 (16)

Let us now consider the problem (11) — optimization of the mismatch of the model criteria. In our case, this problem has the following form:

minimize p according to the set of variables {x1, x2, x3,p} subject to x1 > 0, x2 > 0, x3 > 0, p > 0 and a1x1 + a2x2 + a3x3 < b

In addition, the following inequalities must be satisfied:

(17)

xk > Fk-p Vk = 1,2,3

This problem can also be solved by using the method of smooth penalty functions with the same P(t,s), for which it will be convenient to introduce (in addition to the previously defined) the notation Yk = Fk -p-xk Vk = 1,2,3 . In this case, the auxiliary function (12) will be:

E (t, x,p,u ) = —p — t exp I — I-£t exp I — I-^t exp

k=1

y

j=1

v t /

(18)

The conditions of stationarity for the auxiliary function (13) in (18) take the form of a system of equations:

|p= -1+^ (vl-it exp fv]=0

k=1

dE = expf+ expakexp| ^ | = 0 Vk = 1,2,3

(19)

The solution of the system (19) is denoted by p(u) and x (u), and the function below is used as the smooth approximation of the dependency p** (u):

E (u ) = -E (t, p(u), x (u ),u) (2Q)

The derivatives of this function with respect to the components of the parameter vector u are found from formulas (14) and for the considered example have the form:

(E(u 1 = exP (*). +iexp f ^

= exp

y 4

I 3

da, db

p k=1 3

ix,

v j=1 dup dup j

k=1

t J 5u„

K ] dF

2H f J^ vp = ^

The values of the derivatives of the functions Fk (u) are found from the formulas (16).

Now, for searching for the stationary points of the function E (u) in the space of parameters, it is possible to use any iterative method that improves the directions which can be found by using the derivatives of the first order. As an example, for the solution of problem (17) we apply standard gradient scheme of steepest ascent.

Let t = 0,1,2,... — iteration number. Then this scheme can be described by using the relations:

ut+l = ut + stwt

p p p

where

wp = p N

- ( E (u

p = 1,2.

The value of the norm of the gradient is calculated from the usual formula of the orthonormal basis:

N^, =

(E (')). | +[(*>'))

and the quantity s' — the step along the improving direction for each iteration by the dichotomy method.

The results are shown in tables 1 and 2. Remarks on the use of algorithms of second order

In conclusion, we consider the method of finding elements of the Hessian matrix for the function E (u), the knowledge of which will allow us to use second-order methods in the search for stationary points.

p

Table 1

t u1 u2 E p

0 0.700000000 1.600000000 0.580923855 0.545812501

1 0.700000000 1.200000000 0.633041421 0.596126653

2 0.900876900 1.159092100 0.654535635 0.621144210

3 0.941239800 1.000136600 0.660271356 0.626927551

4 0.981719000 1.018621200 0.661487390 0.628152833

5 1.000234400 1.000052500 0.661620557 0.628286990

6 1.000071000 0.999946090 0.661620583 0.628287016

Table 2

t Ngrad Wi W2 s

0 0.231363725

1 0.131142557

2 0.076755604

3 0.050581430

4 0.010207109

5 2.43444 *10A-4

6 3.73575 *10A-5

-1.8973 * 10A-9 0.979887531 0.246115279 0.909645878 0.706089900 -0.837971176 -0.922621240

-1.000000000 -0.199550560 -0.969240563 0.415384614 -0.708122202 -0.545714493 0.385707203

0.400000000 0.205000000 0.164000000 0.044500000 0.026222500 0.000195000 0.000076500

Applying the rules for differentiating a composite function to the function (20), we obtain:

(Wi \\ d2E " d2E dXj d2E dp ri ,

(E (u)) =-+ V--L+--— p,q = 11,r I

v 'upuq du du j=1 du dXj du du dp du

(21)

The second partial derivatives are calculated directly at the point {x , p}, and the first derivatives, i.e. dXj and — are found according to the well-known implicit function theorem from the system

du

du„

"q - "q

of linear equations:

d2E dp " d2E dxi d2E n

—2—— +y--+-= 0

dp2 duq dpdXi duq dpduq

d2E dp " d2E 8Xj d2E . w. ri , --—+ >--L +-= 0 Vj = 11,« I

dx .dp du ^^ pv ^v ^v PJI ^

q i=1 dXj dXi 5uq dXj5uq

which is obtained by successively differentiating the stationarity conditions (13) according to the variables p and xL j = [1,n].

Finally, it must be mentioned that to calculate the derivatives (21) it is also necessary to know the values of the second derivatives of the functions Fk (u). These values can be found (similar to the

one used above) by the method from formulas (10) and conditions (9) — the stationarity of the functions Ak (t,x, u).

References

Fiacco, A. V., & McCormick, G. P. (1968). Nonlinear Programming: Sequential Unconstrained Minimization, Techniques.

New York, NY: John Wiley and Sons. Lotov, A. V., & Pospelov, I. I. (2008). Mnogokriterial'nye zadachi prinyatiya reshenii [Multi-criterial decision-making

tasks]. Moskva, Rossia: MAKS Press. Umnov, A. E. (1975). Metod shtrafnykh funktsii v zadachakh bol'shoi razmernosti [The method of penalty functions in problems of large dimension]. Zhurnal vychislitel'noi matematiki i matematicheskoi fiziki, 15(6), pp. 1451-1463. Zhadan, V. G. (2014). Metody optimizatsii. Chast' 1. Vvedenie v vypuklyi analiz i teoriyu optimizatsii [Methods of optimization. Part 1. Introduction to convex analysis and optimization theory]. Moskva, Rossia: MFTI.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Оптимизация формы множества Парето в задачах многокритериального программирования Була Аксел Катенди1, Е.А. Умнов2, А. Е. Умнов2

В работе рассматривается схема использования метода гладких штрафных функций для исследования зависимости решений задач многокритериальной оптимизации от параметров. Приводится описание алгоритмов, основанных на методе гладких штрафных функций, решения задачи оптимизации по параметрам уровня согласованности целевых функций и выбора соответствующей формы множества Парето.

Ключевые слова: задача многокритериального параметрического программирования; множество Парето; метод гладких штрафных функций; задача оптимизации по параметрам JEL classification: C61, C65

u Московский физико-технический институт (Государственный университет), Москва, Россия 1 [email protected], [email protected]

i Надоели баннеры? Вы всегда можете отключить рекламу.