<<

стр. 3
(всего 3)

СОДЕРЖАНИЕ

x1 x2 = x3 x4 , b12 > ?12 + ?34 ; x1 x3 = x2 x4 , b13 > ?13 + ?24 ;
x1 x4 = x2 x3 , b14 > ?14 + ?23 . (11.29)
В реальных задачах влияние тройных взаимодействий обычно равно
нулю. Следовательно, генерирующее соотношение (11.26) следует ис-
пользовать, если наибольший интерес представляют оценки для ли-
нейных эффектов.
Для соотношения (11.27) определяющим контрастом будет
I = x1 x3 x4 . (11.30)
Тогда

52
x1 = x3 x4 , b1 > ?1 + ?34 ; x2 = x1 x2 x3 x4 , b2 > ? 2 + ?1234 ;
x3 = x1 x4 , b3 > ?3 + ?14 ; x4 = x1 x3 , b4 > ? 4 + ?13 ;
x1 x2 = x2 x3 x4 , b12 > ?12 + ? 234 ; x2 x3 = x1 x2 x4 , b23 > ? 23 + ?124 ;
x2 x4 = x1 x2 x3 , b24 > ? 24 + ?123 . (11.31)
Следовательно, дробную реплику с генерирующим соотношением
x4 = x1x3 следует использовать, если наибольший интерес представля-
ют эффекты парных взаимодействий.
В общем случае число опытов в дробной реплике должно удовле-
творять следующему соотношению:
k + 1 ? N < 2k , (11.32)
где k — число факторов. Если число опытов равно числу определяе-
мых коэффициентов в линейном уравнении регрессии (N = k + 1),
дробная реплика представляет собой линейный насыщенный план, для
которого все линейные эффекты смешаны с эффектами взаимодейст-
вия. Число степеней свободы остаточной дисперсии в таких планах
равно нулю, поэтому для проверки адекватности линейного уравнения
необходимо проведение дополнительных опытов.
Итак, рассмотренные двухуровневые планы ПФЭ 2k и ДФЭ 2k-1 об-
ладают следующими свойствами: вычисления просты; все коэффици-
енты регрессии определяются независимо друг от друга и с одинако-
вой и минимальной дисперсией; каждый коэффициент рассчитывается
по результатам всех опытов.




53
ЛЕКЦИЯ 12
Оптимизация методом крутого восхождения по поверхности отклика.
Описание функции отклика в области, близкой к экстремуму. Компози-
ционные планы Бокса-Уилсона. Ортогональные планы второго порядка,
расчет коэффицентов уравнения регрессии. Метод последовательного
симплекс-планирования.

12.1. Оптимизация методом крутого восхождения по поверхно-
сти отклика
Задача оптимизации сводится к опытному определению такого со-
четания уровней k факторов (координаты точки в (k+1)-мерном фак-
торном пространстве), при котором достигается максимальное (мини-
мальное) значение выходного параметра y (или нескольких парамет-
ров), т. е. функция отклика системы
y = ? ( x1 , x2 ,..., xk )
принимает экстремальное значение.
Рассмотрим случай, когда на систему оказывают влияние только
два фактора (х1 и х2 в безразмерном масштабе). Построим контурные
сечения y = const поверхности отклика при k = 2 (рис. 3 а).




Рис. 3. Движение по поверхности отклика (а) к экстремуму
в традиционном эксперименте и в методе крутого восхождения (б)

Поиск экстремальной точки поверхности отклика в традиционном
эксперименте проводится следующим образом. В точке L с известным
значением у фиксируется один из факторов, например х1, и начинается
движение из этой точки вдоль оси х2. Движение по х2 продолжается до
тех пор, пока не прекращается прирост y (рис. 3 б). В точке М с наи-
лучшим значением выходного параметра фиксируется фактор х2 и на-
чинается движение в направлении оси х1. В точке N со следующим


54
наилучшим значением y снова фиксируется х1 и начинается движение
по х2 и т. д. Очевидно, что путь к экстремуму по ломаной кривой
LMNR (рис. 3 б) не является оптимальным.
Кратчайшим, наиболее крутым путем достижения экстремума бу-
дет движение из точки L по градиенту перпендикулярно изолиниям
y = const (на рис. 3 б этот путь показан пунктирной линией). Для рас-
сматриваемого случая градиент функции отклика равен
? ?? ?> ? ?? ?>
grad ? = ?
??x ? i +??x ? j , (12.1)
? ? ?
? 1? ? 2?
> >
где i и j — орты координатных осей. Предполагается, что функция
? непрерывна, дифференцируема и не имеет особых точек.
Для реализации метода крутого восхождения Бокс и Уилсон пред-
ложили шаговый метод движения по поверхности отклика. В окрест-
ности точки L ставится эксперимент для локального описания поверх-
ности отклика линейным уравнением регрессии:
?
y = b0 + b1 x1 + b2 x2 . (12.2)
Движение из точки L начинается в направлении градиента линейного
приближения
? ?
?y ?y
= b1; = b2 . (12.3)
? x1 ? x2
Для случая, представленного на рис. 3 б, выборочные коэффици-
енты при линейных членах в окрестности точки L имеют разные зна-
ки: b1 > 0, b2 < 0, поэтому при движении к максимуму функции откли-
ка значение х1 увеличивается, а х2 уменьшается. Движение по гради-
енту линейного приближения продолжается до тех пор, пока не пре-
кращается прирост y. В точке с наибольшим значением y (центр пла-
на) ставится новая серия опытов и определяется новое направление
движения по поверхности отклика. Такой шаговый процесс продол-
жается до достижения области, близкой к экстремуму.
При постановке опытов величина шага должна быть пропорцио-
нальна произведению коэффициента на интервал варьирования: bj ?zj.
Например, при движении из точки L следующий эксперимент ставит-
ся в точке со значениями x1 и х2, отличающимися от начальных на ве-
личины 2b1?z1 и 2b2?z2 соответственно. В общем случае направление
градиента будет зависеть от выбранного интервала варьирования не-
зависимых факторов. При изменении в n раз интервала варьирования
некоторого j-фактора величина шага для него меняется в n2 раз, так

55
как при этом в n раз изменяется и коэффициент регрессии bj. Инвари-
антными к изменению интервала остаются только знаки составляю-
щих градиента. При увеличении числа рассматриваемых факторов бо-
лее двух оптимизация методом крутого восхождения по поверхности
отклика проводится аналогичным способом.

12.2. Описание функции отклика в области, близкой к экстре-
муму. Композиционные планы Бокса-Уилсона
В области, близкой к экстремуму, (или «почти стационарной об-
ласти») функция отклика существенно нелинейна, поэтому для ее аде-
кватного описания необходимо использовать нелинейные полиномы.
В настоящее время для этой цели наиболее широко применяют поли-
номы второго порядка, для получения которых имеются хорошо раз-
работанные планы эксперимента.
Для описания полиномом второго порядка эксперимента, реализо-
ванного для нахождения оптимальных условий процесса, число опы-
тов N в плане должно быть не меньше числа определяемых коэффи-
циентов в уравнении регрессии второго порядка для k факторов
k k k

? ? ?
?
b jj x 2 .
y = b0 + bj x j + buj xu x j + (12.4)
j
j =1 u , j =1 j =1
u? j

Выборочные коэффициенты (12.4) являются оценками соответствую-
щих коэффициентов уравнения теоретической регрессии:
k k k

? ? ? ? jj x 2 .
m y = ?0 + ? jxj + ?uj xu x j + (12.5)
j
j =1 u , j =1 j =1
u?j

В зависимости от числа рассматриваемых факторов число коэф-
фициентов l уравнения регрессии (12.4) определяется по формуле
(k + 1)(k + 2)
k!
2
l = (k + 1) + k + Ck = 2k + 1 + = , (12.6)
k ? 2)!
2!( 2
2
где Ck — количество сочетаний из k факторов по два, равное числу
эффектов парного взаимодействия.
В области, близкой к экстремуму, становятся значимыми эффекты
парного взаимодействия и квадратичные эффекты. Поэтому то, что
адекватное описание результатов эксперимента требует использова-
ния полиномов второго порядка, может служить признаком нахожде-
ния в почти стационарной области. Близость к этой области можно

56
также установить, поставив дополнительно к ПФЭ 2k или ДФЭ 2k-1 се-
рию опытов в центре плана. Среднее значение результатов этих опы-
тов является оценкой для свободного члена уравнения (12.5):
y o > ?0 . (12.7)
Выборочный коэффицент b0, вычисляемый по формуле
N

?x
1
b0 = 0i y i , (12.8)
N i =1

оценивает сумму свободного и квадратичных членов:
k

??
b0 > ? 0 + . (12.9)
jj
j =1

Поэтому, чем больше разность
k

??
o
(b0 ? y ) > , (12.10)
jj
j =1

тем значимее квадратичные эффекты.
Для описания поверхности отклика полиномами второго порядка
независимые факторы в планах должны принимать не менее трех раз-
ных значений. Эксперимент, в котором каждый из k факторов рас-
сматривается на трех уровнях и реализуются все возможные сочета-
ния уровней факторов, является ПФЭ типа 3k. В качестве примера в
табл. 15 представлена матрица планирования ПФЭ 32.
Проведение ПФЭ 3k требует большого числа опытов, намного пре-
вышающего число определяемых коэффициентов l в уравнении (12.4)
уже при k > 2:


2 3 4
k
3k 9 27 81
6 10 15
l




57
Таблица 15

Матрица планирования ПФЭ 32
№ опыта х1 х2 y
1 –1 –1 y1
2 0 –1 y2
3 +1 –1 y3
4 –1 0 y4
5 0 0 y5
6 +1 0 y6
7 –1 +1 y7
8 0 +1 y8
9 +1 +1 y9

Сократить общее число опытов при условии получения несмешан-
ных оценок для линейных эффектов и эффектов взаимодействия мож-
но с помощью композиционных планов Бокса-Уилсона. Ядро таких
планов при k < 5 составляет ПФЭ 2k, и полуреплика от него при k ? 5.
Если линейное уравнение регрессии оказалось неадекватным экс-
перименту, необходимо:
1) добавить 2k звездных точек, расположенных на координатных
осях факторного пространства. Координаты звездных точек в
общем случае равны
(± ?, 0, ..., 0), ( 0, ± ?, ..., 0), ..., ( 0, ..., 0, ± ?) ,
где ? — расстояние от центра плана до звездной точки, или
звездное плечо;
2) увеличить число экспериментов в центре плана n0.
Общее число опытов в матрице композиционного плана при k ? 4
составляет
N = 2 k + 2k + n0 . (12.11)
Рассмотрим построение композиционных планов на примере k = 2
(рис. 4). Точки 1, 2, 3, 4 образуют ПФЭ 22, точки 5, 6, 7, 8 являются
звездными точками с координатами (± ?, 0) и (0, ± ? ) , координаты n0
опытов в центре плана нулевые — (0, 0).
Композиционный план второго порядка для двух факторов пред-
ставлен в табл. 16, при этом в центре плана выполнена серия из трех
опытов (№ 9–11).




58
Рис. 4. Композиционный план второго
порядка для двух факторов
Таблица 16

Композиционный план второго порядка для двух факторов
x0 x1 x2 x1 x2 x12 2
№ опыта y
x2
1 +1 +1 +1 +1 +1 +1 y1
2 +1 +1 –1 –1 +1 +1 y2
3 +1 –1 –1 +1 +1 +1 y3
4 +1 –1 +1 –1 +1 +1 y4
?2
5 +1 0 0 0 y5
+?
?2
6 +1 0 0 0 y6
–?
?2
7 +1 0 0 0 y7
+?
?2
8 +1 0 0 0 y8
–?
9 +1 0 0 0 0 0 y9
10 +1 0 0 0 0 0 y10
11 +1 0 0 0 0 0 y11

12.3. Ортогональные планы второго порядка, расчет коэффи-
центов уравнения регрессии
Выбор звездного плеча в композиционных планах Бокса–Уилсона
может быть произвольным, однако расчеты коэффициентов уравнения
регрессии при k < 5 существенно упрощаются, если величина плеча
определяется исходя из следующего уравнения:
? 4 + 2k? 2 ? 2 k ?1 (k + 0.5n0 ) = 0 . (12.12)
Значения ?2, определенные по (12.12), приведены в табл. 17.

59
Таблица 17


Значения ?2 для k факторов и n0 опытов в центре плана
k k
n0 n0
2 3 4 2 3 4
1 1.00 1.476 2.00 6 1.742 2.325 2.950
2 1.160 1.650 2.164 7 1.873 2.481 3.140
3 1.317 1.831 2.390 8 2.00 2.633 3.310
4 1.475 2.00 2.580 9 2.113 2.782 3.490
5 1.606 2.164 2.770 10 2.243 2.928 3.66

Выбрав ?, проведем следующее линейное преобразование квадра-
тичных столбцов:
N

?
1
x 'j = x 2 ? x 2 = x 2 ? x2 . (12.13)
j j j ji
N i =1

Композиционные планы, полученные таким образом, называются ор-
тогональными планами второго порядка.
Ортогональный план второго порядка при k = 2 и n0 = 1 представ-
лен в табл. 18. За его основу взят композиционный план для двух фак-
торов (табл. 16) с общим числом опытов N = 9. Величину звездного
плеча определим по табл. 17: ?2 = 1, ? = 1. Средние значения элемен-
тов квадратичных столбцов в табл. 16 равны
9 9

? ?
1 1 1 2
2 2 2
x2i = (4 + 2a2 ) = .
2
x1 = x1i = x2 = (12.14)
9 9 9 3
i =1 i =1

В математической статистике доказывается, что для ортогональных
планов второго порядка все коэффициенты уравнения регрессии оп-
ределяются независимо друг от друга по формуле
N N

? ? x2 ,
bj = x ji yi (12.15)
ji
i =1 i =1

а дисперсии коэффициентов равны
N

?
s 2 (b j ) = sвоспр
2
x2 . (12.16)
ji
i =1




60
Таблица 18


Ортогональный план второго порядка для двух факторов

x1' '
№ опыта y
x0 x1 x2 x1 x2 x2
1 +1 +1 +1 +1 +1/3 +1/3 y1
2 +1 +1 –1 –1 +1/3 +1/3 y2
3 +1 –1 –1 +1 +1/3 +1/3 y3
4 +1 –1 +1 –1 +1/3 +1/3 y4
5 +1 0 0 +1/3 –2/3 y5
+?
6 +1 0 0 +1/3 –2/3 y6
–?
7 +1 0 0 –2/3 +1/3 y7
+?
8 +1 0 0 –2/3 +1/3 y8
–?
9 +1 0 0 0 –2/3 –2/3 y9

Для определения дисперсии воспроизводимости необходимо вы-
полнить серию опытов в центре плана. В результате расчетов по мат-
рице с преобразованными столбцами для квадратичных эффектов
(табл. 18) получаем следующее уравнение:
?
' 2 2 2 2
y = b0 + b1 x1 + b2 x2 + b12 x1 x2 + b11 ( x1 ? x1 ) + b22 ( x2 ? x2 ) . (12.17)
Чтобы перейти к обычной записи уравнения регрессии в виде
?
2 2
y = b0 + b1 x1 + b2 x 2 + b12 x1 x 2 + b11 x1 + b22 x 2 , (12.18)
определим b0 по формуле
' 2 2
b0 = b0 ? b11 x1 ? b22 x2 (12.19)
с дисперсией, равной

() ()
222 222
2 2 '
s (b0 ) = s (b0 ) + x1 s (b11 ) + x 2 s (b22 ) . (12.20)
Значимость коэффициентов проверяется по критерию Стъюдента.
Если
bj
tj = > tp( f ), (12.21)
s (b j )
где tp(f ) — табличное значение критерия Стъюдента для р = 0,05 и
числа степеней свободы дисперсии воспроизводимости, то коэффици-
ент значимо отличается от нуля.



61
Коэффициенты уравнения регрессии, получаемые при помощи ор-
тогональных планов второго порядка, определяются с разной точно-
стью. В случае, когда k ? 4, согласно (12.16) имеем
2
sвоспр
2 '
=
s (b0 ) ;
N
2
sвоспр
s 2 (b j ) = , j = 1, 2, ..., k ;
k 2
2 + 2?
2
sвоспр
s 2 (buj ) = , u , j = 1, 2, ..., k , u ? j ; (12.22)
k
2
2
sвоспр
s 2 (b jj ) = , j = 1, 2, ..., k .
k
x 2 )2 2
x 2 )2 ? 2)( x 2 ) 2
2 (1 ? + 2(? ? + (n0 + 2k
j j j

После исключения незначимых коэффициентов проводится про-
верка адекватности уравнения по критерию Фишера. Уравнение адек-
ватно эксперименту, если
2
sад
F= ? F1? p ( f ад , f воспр ) , (12.23)
2
sвоспр
где F1-p (fад, fвоспр) — критерий Фишера для р = 0,05; fад = N – l — число
степеней свободы дисперсии адекватности (l — число значимых ко-
эффициентов в уравнении регрессии); fвоспр — число степеней свободы
дисперсии воспроизводимости.

12.4. Метод последовательного симплекс-планирования
В рассмотренных выше планах ПФЭ 22 и 23 экспериментальные
точки располагались в вершинах квадрата и куба соответственно. В
качестве экспериментального плана можно также использовать регу-
лярный симплекс. Симплексом в k-мерном пространстве называют вы-
пуклый многогранник, имеющий ровно (k + 1) вершину, каждая из ко-
торых определяется пересечением k гиперплоскостей данного про-
странства. Симплекс называется регулярным, если расстояния между
всеми его вершинами равны. Примерами регулярных симплексов яв-
ляются правильный треугольник в двумерном пространстве и тетра-
эдр в трехмерном.
На практике планирование эксперимента с использованием регу-
лярных симплексов применяется для решения задач оптимизации при



62
движении к почти стационарной области. Для получения регулярного
симплекса проводится линейное преобразование уровней факторов
z j ? zo
j
xj = , (12.24)
?zj

где z o — j-я координата центра плана; ?z j — интервал варьирования
j
по j-фактору.
Оптимизация методом последовательного симплекс-планирования
проводится следующим образом: исходная серия опытов планируется
так, чтобы экспериментальные точки образовывали регулярный сим-
плекс в факторном пространстве. После проведения опытов определя-
ется вершина симплекса, соответствующая наихудшим результатам.
Далее строится новый симплекс, для чего наихудшая точка исходного
симплекса заменяется новой, расположенной симметрично относи-
тельно центра грани симплекса, находящейся против наихудшей точ-
ки. Новая точка вместе с оставшимися точками образует новый сим-
плекс, центр тяжести которого смещен в сторону повышения качества
процесса. После реализации опыта в дополнительной точке опять
проводится выявление наихудшей вершины симплекса и т. д. При
достижении области оптимума, симплекс начинает вращение вокруг
вершины с максимальным значением отклика.
На рис. 5 показаны схемы достижения экстремума поверхности
отклика методами крутого восхождения и симплекс-планирования на
примере зависимости целевой функции y от двух факторов. При оп-
тимизации методом крутого восхождения (рис. 5 а) в окрестности
точки М поставлен ПФЭ 22, движение по градиенту линейного при-
ближения осуществлялось в опытах 5–9. Далее был поставлен новый
ПФЭ 22 (точки 10–13) с центром в точке 7, в которой было получено
наилучшее значение y. Движение по новому градиенту (точки 14–15)
приводит к экстремуму.
При оптимизации методом симплекс-планирования (рис. 5 б) в ис-
ходном симплексе (точки 1–3) худшей точкой оказалась точка 2. Ее
зеркальным отражением относительно с1 — центра грани 1–3 — явля-
ется точка 4. В новом симплексе 1, 3, 4 худшей оказалась точка 1, в
результате ее зеркального отражения получен симплекс 3, 4, 5 и т. д.
Область оптимума достигается при реализации симплекса 9, 10, 11.
Хотя оба рассмотренных метода требуют проведения примерно оди-
накового числа опытов, симплекс-планирование имеет ряд важных




63
Рис. 5. Достижение экстремума поверхности отклика
методами крутого восхождения (а) и симплекс-планирования (б)

преимуществ: при использовании этого метода параметр оптимизации
y может измеряться приближенно, достаточно иметь возможность
проранжировать его величину; можно одновременно учитывать не-
сколько параметров оптимизации; метод не предъявляет жестких тре-
бований к локальной аппроксимации поверхности отклика уравнени-
ем регрессии.
На практике рекомендуется ориентировать исходный симплекс в
факторном пространстве следующим образом: центр симплекса сов-
падает с началом координат, одна из вершин лежит на координатной
оси, а остальные располагаются симметрично относительно коорди-
натных осей, плоскостей и гиперплоскостей (в многомерном случае).
Тогда координаты вершин симплекса при k = 5 задаются матрицей
? x1 x5 ?
x2 x3 x4
?? x x5 ?
x2 x3 x4
?1 ?
? 0 ? 2 x2 x5 ?
x3 x4
X =? ?. (12.25)
0 ? 3 x3
0 x4 x5 ?
?
?0 x5 ?
0 ? 4 x4
0
? ?
0 ? 5 x5 ?
0 0 0
?
При k < 5 координаты вершин симплекса определяются частью мат-
рицы (12.25), при этом число столбцов равно числу факторов, а число
строк равно k + 1; при k > 5 для каждого добавленного фактора в мат-
рицу (12.25) добавляются соответствующие столбец и строка. В об-

64
щем случае число опытов в симплексной матрице для k независимых
факторов N = (k + 1) равно числу коэффициентов линейного уравне-
ния регрессии, т. е. симплексные планы являются насыщенными.
Если длину стороны симплекса принять равной 1, то
1
xj = . (12.26)
2 j ( j + 1)
Для практического использования матрицы (12.25) ее числовые эле-
менты заранее подсчитаны по формуле (12.26)
? 0.5 0.129?
0.289 0.204 0.158
?? 0.5 0.129?
0.289 0.204 0.158
? ?
? 0 ? 0.578 0.129?
0.204 0.158
X =? ?. (12.27)
0 ? 0.612
0 0.158 0.129?
?
?0 0.129?
0 ? 0.632
0
? ?
0 ? 0.645?
0 0 0
?
План эксперимента в безразмерном масштабе для k факторов со-
стоит из k столбцов и k + 1 строки матрицы (12.27). Коэффициенты
уравнения линейной регрессии
k
?
?b x
y = b0 + (12.28)
j j
j =1

вычисляются следующим образом:
N N

?y , ?x
1
b0 = bj = 2 ji yi . (12.29)
i
N i =1 i =1

Если в одной из вершин симплекса поставить серию параллельных
опытов и рассчитать дисперсию воспроизводимости, то выборочные
дисперсии коэффициентов определяются по формуле
2
sвоспр
s 2 (b j ) = 2
= 2 sвоспр . (12.30)
N

?x 2
ji
i =1

Следует отметить, что коэффициенты уравнения регрессии, полу-
ченные по симплексному плану, определяются с меньшей точностью
по сравнению с коэффицентами, полученными при реализации ПФЭ
2k и ДФЭ 2k-1, для которых


65
s 2 (b j ) = sвоспр / N .
2
(12.31)
После реализации исходного симплекса требуется провести отра-
жение наихудшей точки относительно центра противоположной гра-
ни. Координаты отраженной точки равны:
x (jk + 2) = 2 x (jc ) ? x (jl ) , j = 1, 2, ..., k , (12.32)

где x (l ) — j-я координата наихудшей точки; x (jk + 2) — j-я координата
j

новой точки, получаемой в результате отражения; x (c ) — j-я коорди-
j
ната центра противоположной грани, определяемая по формуле
k +1

?
1
x (jc ) x (ji ) , i ? l ,
= (12.33)
k
i =1

где x (i ) — j-я координата i-й вершины симплекса (i = 1, 2, …, k+1).
j
Координаты центра оптимального симплекса (точка S) в почти
стационарной области находятся следующим образом:
k +1

?
1
x (jS ) = x (ji ) . (12.34)
(k + 1)
i =1




66
Учебное издание

Блохин Андрей Викторович

ТЕОРИЯ
ЭКСПЕРИМЕНТА

КУРС ЛЕКЦИЙ
В двух частях
Часть 2
В авторской редакции

Художник обложки А. А. Федорченки
Технический редактор Г. М. Романчук
Корректор Р. П. Кадырко

Ответственный за выпуск А. В. Блохин


Подписано в печать 20.11.2002. Формат 60x84/16. Бумага офсетная.
Гарнитура Тайме. Печать офсетная Уел печ. л. 3,95. Уч.-изд. л. 3,62.
Тираж 100 экз. Зак.58

Белорусский государственный университет.
Лицензия ЛВ №315 от 14.07.98.
220050, Минск, проспект Франциска Скорины, 4.

Отпечатано с оригинала-макета заказчика.
Республиканское унитарное предприятие
«Издательский центр Белорусского государственного университета».
Лицензия ЛП№ 461 от 14.08.2001.
220030, Минск, ул. Красноармейская, 6.

<<

стр. 3
(всего 3)

СОДЕРЖАНИЕ