Обработал данные третьего плана многофакторного планирования и получил опять результаты, которые позволяют сделать только один вывод скорость распространения гравитации должна быть больше чем 200 скоростей света. Это меня заставило задуматься о том все ли у меня в порядке с теоретическим обоснованием этого моего исследования, т.к., проводя ранее подобные исследования, я обычно за два шага уверенно приходил в область оптимума, а здесь сделал три шага и оптимума не видно. И я стал даже задумываться о том, что может быть методы многофакторного планирования по тому критерию оптимизации, что я использую, не очень подходят для этого. Ведь до этого исследования, если не считать случая по оптимизации коэффициентов в формуле Планка, я оптимизировал параметры систем по отклику системы, а не по разнице между откликом и заданным оптимальным значением (в нашем случае в формуле (3), которую я приводил выше, между расчетными, т.е. полученными при вычислительном эксперименте на модели, YRas(I, J, U) и наблюдаемыми YNab(I, J) значениями вековых смещений в U-ом эксперименте для I ой планеты и J-го параметра). В дальнейшем этот критерий я буду называть dY в противовес критерию Y, который обычно применяется при многофакторном планировании, и где, при проведение натурных экспериментов, Yu0(U)= YNab(U), а, при проведение вычислительных экспериментов, Yu0(U)= YRas(U).

Yu0(U) = SUMi,j ( kVesa(I, J) * Abs ((YRas(I, J, U) - YNab(I, J)) / YNab(I, J)) / 100) (3)

Вообще то в книге С.В.Мельников, В.Р.Алешкин, П.М.Рощин Планирование эксперимента в исследованиях сельскохозяйственных процессов Л. Колос 1980 на стр. 45 для расчета комплексного критерия приводится формула подобная моей формуле (3). Там только используется не абсолютное значение относительной разности между откликом системы и оптимальным значением, а квадрат этой разности, т.е. критерий dY^2, но, как пишут авторы, это делается только для того, чтобы разность была всегда положительной. Я же в своей формуле использовал абсолютное значение этой разницы, т.е. принципиальных отличий от их формулы у меня нет и, следовательно, я могу смело использовать свой критерий оптимизации dY для оптимизации параметров Солнечной системы. Но одно дело, что у них там написано, а другое дело то, что я вижу. Да и мой собственный опыт с коэффициентами в формуле Планка (1t) не очень показательный, т.к. оптимизировал я там только 3 коэффициента, а 4-ый фактор (температура излучения) мною принудительно задавался для повышения качества информации при проведение вычислительных экспериментов по почти D-оптимальному плану Бокса для четырех факторов. К тому же второй коэффициент (показатель степени при частоте излучения v) мне надо было не столько оптимизировать, сколько подтвердить, что он равен 3, как это следовало из формулы Вина. Таким образом, я по большому счету оптимизировал только 2 параметра (фактора) и, следовательно, у меня могли быть только двойные смешанные взаимодействия, а это прекрасно воспроизводится полиномом 2-ой степени (2), который я получаю после обработки данных вычислительных экспериментов.

А кто может ответить на вопрос есть ли в нашей системе, которую мы исследуем, смешанные взаимодействия выше парных, т.е. тройные или четверные. А может быть даже есть и не только линейные взаимодействия, но и квадратичные. К сожалению, ответить на эти вопросы никто не может. Да, наверное, никто не сможет ответить и на то, как это скажется на описание поверхности отклика при таких условиях. По этому, я на всякий случай (не очень доверяя всему, что написано в учебниках) решил провести маленькое исследование по оптимизации по критериям dY и dY^2 параметров простейших математических выражений, которые будут имитировать поведение различных систем. И первым делом я решил взять чуть ли не самый сложный случай с четверным взаимодействием, где вдобавок одно взаимодействие еще и не линейно, т.е. всем Вам известный закон тяготения Ньютона (3t) и попробовать оптимизировать его параметры. В принципе, мы можем с законом тяготения провести и натурные эксперименты. Правда, не с самим законом тяготения для масс, а с законом тяготения для зарядов (закон Кулона), где даже аналог гравитационной постоянной можем изменять, распологая различные диэлектрики между зарядами. Но речь сейчас идет не о том, можем ли мы воспроизвести эксперименты на реальном объекте или на его модели, а о том, можем ли мы, уже даже зная аналитическую формулу, отражающую отклик системы на наши воздействия на нее, чисто с математической точки зрения получить оптимальные значения системы по примененному мною критерию dY, т.е. по разнице между откликом системы и известным оптимальным значением.

Может возникнуть вопрос а зачем вообще надо проводить исследования для получения аппроксимации (2), если у нас уже есть аналитическая формула закона тяготения. А затем, что, мы сейчас просто проверяем на что способны методы многофакторного планирования, чтобы заранее знать, что от них ожидать. Ведь когда мы исследуем какую то сложную систему, то нам надо проводить натурные или вычислительные эксперименты, чтобы получить хотя бы уравнение регрессии (2), т.к. никакие аналитические выражения для критерия оптимизации при исследование самого объекта нам не известны вообще, а аналитическая формула, по которой вычисляется критерий оптимизации в моделях объекта, даже если и удастся такую получить в развернутом виде, может уместиться только на десятках или сотнях страниц, что делает ее не пригодной для аналитических методов оптимизации. А уравнение регрессии (2), т.е. полином 2-го порядка, который мы получаем при многофакторном планировании, очень удобен для этого и по этому мы и постараемся его получить по критерию dY для тестируемых систем. А т.к. в программе Solsys5 у меня по формуле (3) рассчитывается значение комплексного критерия оптимизации (целевой функции) в каждом из 24 экспериментов, а отклик системы в наших тестовых примерах определяется не по комплексному критерию, то мы можем, для оптимизации параметров по критерию dY в тестовых выражениях (3t10t), формулу (3) упростить до выражения (4)

Yu0(U) = Abs((YRas(U) - Yopt) / Yopt) (4)
Где: Yu0 (U) относительная разница между расчетным YRas(U) и оптимальным Yopt значением отклика системы, поведение которой имитирует одна из формул (3t10t), в U-ом эксперименте.

Результаты оптимизации параметров в формуле тяготения, по примененному мною критерию оптимизации dY, получились удручающие, т.к. аппроксимация критерия оптимизации, полученным уравнением регрессии (2), не лезла ни в какие ворота. Да Вы сами взгляните на полученные значения критерия оптимизации dY с использованием формулы закона тяготения (3t) и эти же значения по полученному уравнению регрессии (2) на нижеприведенном рисунке (верхняя часть рисунка), где маленькие синие кружки это критерий оптимизации рассчитанный с использованием формулы (3t) для выражения (3) в 24 экспериментах плана Бокса (номер соответствует абсциссе), а большие синие кружки это критерий оптимизации, рассчитанный по уравнению регрессии (2), для тех же значений параметров, что и в соответствующем эксперименте. При этом все факторы X1 X4, при выполнение плана Бокса, на нулевом уровне были равны единице, а интервалы их варьирования были 0,5, а значение Yopt бралось равным 1, т.е. я принимал, что оптимальные значения параметров X1 X4 в формуле (3t) равны 1 и получалось, что Yopt = 1*1*1/1=1. Для меня такой результат аппроксимации был большой неожиданностью, т.к. ранее в своей книге я сам закон тяготения, т.е. по критерию Y, аппроксимировал уравнением (2) и никаких проблем по качеству аппроксимации тогда не было. Я тут же аппроксимировал полученные значения Yu0(U) полиномом (2) и по критерию Y и выяснилось, что на этот раз у меня получились значительные погрешности в аппроксимации (синие маленькие и большие кружки на нижней части рисунка).

Yu = k0 + k1*X1 + k2*X2 + k3*X3 + k4*X4 +
+ k5*X1*X2 + k6*X1*X3 + k7*X1*X4 + k8*X2*X3 + k9*X2*X4 + k10*X3*X4 +
+ k11*X1^2 + k12*X2^2 + k13*X3^2 + k14*X4^2 (2)
где Yu (или dYu или dYu^2) критерий оптимизации, который надо минимизировать, X1 X4 - оптимизируемые параметры, а k0 k14 коэффициенты, которые мы получаем методом наименьших квадратов при статистической обработке значений Yu0(U) полученных в 24 экспериментах при разных значениях параметров X1 X4.

http://ser.t-k.ru/Ris/3t_Y_dY.gif (зеркало http://modsys.narod.ru/Ris/3t_Y_dY.gif)

Стал разбираться в чем дело и выяснил, что когда я ранее аппроксимировал закон притяжения, то задавал интервалы варьирования параметров примерно 20%, а сейчас задал их 50%. Уменьшил интервал варьирования до 20% и на сей раз опять получил приличный результат (черные кружки), т.е. получается, что все дело не в самом критерии оптимизации, а в интервалах варьирования, тем более, что при интервалах варьирования 80% (зеленые кружки), результат получился еще хуже, чем при 50%. И по критерию dY (в верхней части рисунка) тоже получается, что, чем больше интервал варьирования, тем хуже аппроксимация экспериментальных данных, но здесь уже не все так однозначно, т.к. получается, что и при интервале варьирования 20% аппроксимация экспериментальных данных получается не на много лучше. По этому я решил продолжить исследование с использованием других математических выражений имитирующих поведение системы.

YRas(U) = X1*v^X2* exp(-X3*v/X4) (1t)
YRas(U) = u (2t)
YRas(U) = X1 * X2 * X3 / X4^2 (3t)
YRas(U) = X1^2 * X2^2 (4t)
YRas(U) = X1 * X2 * X3 * X4 (5t)
YRas(U) = X1 + X1^2 (6t)
YRas(U) = X1 + X2 + X3 + X4 (7t)
YRas(U) = X1*X2 + X1*X3 + X1*X4 + X2*X3 + X2*X4 + X3*X4 (8t)
YRas(U) = X1^2 + X2^2 + X3^2 + X4^2 (9t)
YRas(U) = X1 + X2 + X3 + X4 + X1^2 + X2^2 + X3^2 + X4^2 +
+ X1*X2 + X1*X3 + X1*X4 + X2*X3 + X2*X4 + X3*X4 (10t)

При проведение вычислительных экспериментов по всем этим выражениям я принимал, что оптимальные значения всех параметров в этих формулах равны единице и находил сначала оптимальное значение отклика системы Yopt, а затем задавал в соответствие с планом различные значения параметров X1 X4 и, вычислив значение YRas(U), находил значение критерия оптимизации в U-ом эксперименте. Значение критерия Y я брал равным YRas(U), критерия dY я вычислял по формуле (4), а критерий dY^2 определял возведя критерий dY в квадрат. Затем, обработав данные по 24 значениям этих критериев Yu0(U), я получал коэффициенты k0 k14 для аппроксимации (2) по которой вычислял значение Yu для тех же значений параметров X1 X4, что были в 24 экспериментах по плану эксперимента, и сравнивал их с 24 значениями Yu0(U). А по результатам сравнения я выставлял оценку почти D-оптимальному плану Бокса по качеству аппроксимации экспериментальных данных выражением (2) по различным критериям оптимизации. Сравнение я проводил графически определяя попало ли значение Yu0(U), которое на приведенных выше рисунках отражено в 24 экспериментах маленькими кружками, внутрь большого кружка отражающего значение Yu. При этом графический масштаб для вывода данных выбирался так, чтобы все данные от минимального до максимального значения по ординате укладывались в интервале от 1 до 7 сантиметров. А оценки качеству аппроксимации значений критериев оптимизации Y, dY и dY^2, полученных по плану Бокса полиномом (2), я производил по пятибальной шкале, определяя наилучшую оценку по наихудшим результатам, а затем, полученные результаты, оформил в виде таблицы. Критерии оценок были такими

5 баллов все 24 маленьких кружка, т.е. значения Yu0(U), находятся внутри соответствующего большого кружка, центр которого соответствует ординате значения Yu. При этом диаметр большого кружка в два раза больше диаметра маленького.
4 балла - все 24 маленьких кружка или находятся внутри больших кружков или хотя бы касаются его с наружной стороны. При этом я также указываю в таблице, рядом с оценкой, в знаменателе количество экспериментов, по которым была выставлена эта наихудшая оценка. Так, если в 23 случаях маленькие кружки находятся внутри больших, а в одном случае маленький кружок пересекается с большим или касается его с наружной стороны, то будет указана оценка 4/1.
3 балла если все маленькие кружки находятся хотя бы на расстояние 1-го диаметра большого кружка от его окружности. При этом, если 20 маленьких кружков находятся внутри больших, 1 кружок пересекается с большой окружностью, а три маленьких кружка находятся на расстояние 1-го большого диаметра от его окружности, то оценка будет 3/3. Таким образом, в знаменателе указывается только количество кружков не прошедших по более высоким оценкам.
2 балла если все маленькие кружки находятся хотя бы на расстояние 2-х диаметров большого кружка от его окружности.
1 балл - если хотя бы один маленький кружок находится на расстояние более 2-х диаметров большого кружка от его окружности.

Таблица 15. Оценки качества аппроксимации полиномом 2-ой степени различных критериев оптимизации при проведение вычислительных экспериментов по почти D-оптимальному плану Бокса, когда эти критери расчитывались по разным математическим выражениям (имитирующим поведение системы) и при разных интервалах варьирования оптимизируемых факторов (параметров), которые на нулевых уровнях принимали значение 1.

Формула_________критерий_Y________критерий_dY________критерий_dY^2
Интервалы_+/-0,2_+/-0,5_+/-0,8_+/-0,2_+/-0,5_+/-0,8_+/-0,2_+/-0,5_+/- 0,8
3t__________4/2___2/1___1/2______3/3___2/2___1/2_____1/2___1/3___1/6
4t________ ___5____4/4___3/4______3/4___3/4___3/4_____2/4___2/4___2/4
5t___________5___3/10___1/1__ ____2/1___1/2___1/6_____1/2___1/6___1/6
6t___________5____5_____5________5_____5_____5__ ____5_____5_____5
7t___________5____5_____5______3/16___3/16___3/16____5_____5_____5
8t___________5____5_____5________2/1___2/1___2/1_____3/2___2/2___1/2
9t___________5___ _5_____5_______3/11___3/5___3/6_____4/8___4/9___3/6
10t__________5____5_____5______3/16_ __3/15___3/15____4/2___3/2___2/2

Что можно сказать о полученных данных. То, что в тестах 6t10t по критерию Y результаты будут отличными я и не сомневался, т.к. для аппроксимации именно таких взаимодействий полиномы 2-го порядка и приспособлены, но вот почему по критериям dY и dY^2 получились такие плачевные результаты это для меня не понятно. Ясно было и то, что в тестах 3t5t результаты будут хуже чем в тестах 6t10t, т.к. линейные взаимодействия выше парных в примененном нами плане Бокса смешаны с другими взаимодействиями и по этому их нельзя выделить (про парное квадратичное взаимодействие затрудняюсь сказать что то определенное), но и в этом случае я никак не ожидал, что по критериям dY и dY^2 будут такие плохие результаты. По этому остается только надеяться на то, что в исследуемой нами Солнечной системе тройных и выше взаимодействий не будет, а то нам придется очень долго мучиться пока мы найдем область оптимума. И я уже подумываю о включение в 6-ю версию программы ротатабельного плана, у которого информация более равномерно размазана по всем направлениям, т.е. одинаковая на всех направлениях на равных расстояниях от центра плана, а, учитывая то, что звездное плечо у этого плана равно 2, это позволит прощупать факторное пространство на ту же глубину при меньшем значении интервала варьирования для центральной части плана, и по этому этот план может оказаться нам очень полезен именно при поиске области оптимума.

Однако, давайте наконец то посмотрим насколько удачно наш полином (2) аппроксимирует смоделированные нами процессы в Солнечной системе по критериям dY и dY^2. На нижеприведенном рисунке (в верхней части) показаны значения Yu0(U) полученные после обработки данных вычислительных экспериментов на классической математической модели Солнечной системы с учетом скорости распространения гравитации по критерию dY по формуле (3) (маленькие синие кружки) и по критерию dY^2 (маленькие красные кружки) и соответствующие им значения Yu (большие кружки) полученные по уравнениям (2). Причем в уравнение (3) значения рассчитаны по комплексному критерию, т.е. с учетом вековых смещений всех J-х параметров для Меркурия, т.е. J=1 4, а I=1 и с весовыми коэффициентами kVesa(I, J)=100. Как видим по критерию dY аппроксимацию можно оценить на твердую тройку, а если бы не 23 и 24 эксперименты, где у нас был очень большой интервал варьирования скорости гравитации, то возможно бы получилась и четверка. Кстати и интервал варьирования скорости VZsys тоже великоват, т.к. в экспериментах 17 и 18, 19 и 20, 21 и 22, 23 и 24 разность значений критерия между этими парными экспериментами должна быть примерно равна. И по критерию dY^2 аппроксимацию тоже можно оценить на тройку, а если уменьшить интервалы варьирования VZsys и Vgr то тоже может получиться хорошая оценка.

http://ser.t-k.ru/Ris/Plan3.gif (зеркало http://modsys.narod.ru/Ris/Plan3.gif)

Таким образом, быстрее всего, в моделируемой нами системе нет взаимодействий выше парных и аппроксимацию в виде полинома 2-ой степени можно использовать для поиска оптимума. Кстати, если кому-то показалось, что методы многофакторного планирования не очень надежны, то хочу Вас разочаровать, т.к. в основе этих методов лежит очень мощный математический аппарат, хотя и зародились эти методы для исследований в сельском хозяйстве, т.е. в области не очень блещущей математическими изысками. А в доказательство хочу Вам показать маленький фокус, на котором я всегда тестирую все свои новые программы, где применяю многофакторное планирование. На нижней части второго рисунка Вы видите значения Yu0(U) по критерию Y для теста 2t и эти же значения полученные по уравнению регрессии (2) для Yu. А фокус заключается в том, что не зависимо от уровней варьирования факторов мы принимаем значение Yu0(U) в очередном эксперименте просто равным номеру этого эксперимента, т.е. получается, что уравнением (2) мы аппроксимировали почти полный хаос, но результат то получился вполне приличным.

Я конечно же понимаю, что некоторые болельщики за ту или иную теорию уже заждались, когда же я объявлю окончательный результат матча, а я тут рассказываю про скучные двойные взаимодействия и звездные плечи, а им гораздо интереснее было бы услышать про двойные звезды. Но, к сожалению, наука это не сказка и быстро только эти самые сказки сказываются, но не скоро дело делается, а т.к. вопрос мы рассматриваем очень сложный то до окончания этого матча пока еще очень далеко. Хотя, не скрою, что после выполнения третьего плана, я начал уже склоняться к тому, что классическая модель с учетом скорости распространения гравитации (если не в этом матче, то в этом эпизоде уж точно) проиграла. Правда, это, конечно же, не означает, что в этом эпизоде победила ОТОшная модель, т.к. я ее в этих условиях вообще еще не тестировал. А вот, учитывая проведенное мною исследование по описанию систем имитируемых формулами 3t10t, можно пока еще все трудности по нахождению оптимума объяснить именно трудностью описания поверхности отклика системы по критериям dY или dY^2. И хотя, конечно же, как многие догадались, я болею в этом матче именно за классическую модель, но я здесь в лучшем случае арбитр и по этому кто победит покажет только игра. А я сейчас сделаю описание к 5-ой версии программы Solsys и выложу ее для скачивания, чтобы самые нетерпеливые болельщики сами смогли провести на ней вычислительные эксперименты и попытаться найти оптимум для наших параметров системы (если, конечно же, он там есть для классической модели Солнечной системы).

С наилучшими пожеланиями Сергей Юдин.

Обработал данные третьего плана многофакторного планирования и получил опять результаты, которые позволяют сделать только один вывод скорость распространения гравитации должна быть больше чем 200 скоростей света. Это меня заставило задуматься о том все ли у меня в порядке с теоретическим обоснованием этого моего исследования, т.к., проводя ранее подобные исследования, я обычно за два шага уверенно приходил в область оптимума, а здесь сделал три шага и оптимума не видно. И я стал даже задумываться о том, что может быть методы многофакторного планирования по тому критерию оптимизации, что я использую, не очень подходят для этого. Ведь до этого исследования, если не считать случая по оптимизации коэффициентов в формуле Планка, я оптимизировал параметры систем по отклику системы, а не по разнице между откликом и заданным оптимальным значением (в нашем случае в формуле (3), которую я приводил выше, между расчетными, т.е. полученными при вычислительном эксперименте на модели, YRas(I, J, U) и наблюдаемыми YNab(I, J) значениями вековых смещений в U-ом эксперименте для I ой планеты и J-го параметра). В дальнейшем этот критерий я буду называть dY в противовес критерию Y, который обычно применяется при многофакторном планировании, и где, при проведение натурных экспериментов, Yu0(U)= YNab(U), а, при проведение вычислительных экспериментов, Yu0(U)= YRas(U).

Yu0(U) = SUMi,j ( kVesa(I, J) * Abs ((YRas(I, J, U) - YNab(I, J)) / YNab(I, J)) / 100) (3)

Вообще то в книге С.В.Мельников, В.Р.Алешкин, П.М.Рощин Планирование эксперимента в исследованиях сельскохозяйственных процессов Л. Колос 1980 на стр. 45 для расчета комплексного критерия приводится формула подобная моей формуле (3). Там только используется не абсолютное значение относительной разности между откликом системы и оптимальным значением, а квадрат этой разности, т.е. критерий dY^2, но, как пишут авторы, это делается только для того, чтобы разность была всегда положительной. Я же в своей формуле использовал абсолютное значение этой разницы, т.е. принципиальных отличий от их формулы у меня нет и, следовательно, я могу смело использовать свой критерий оптимизации dY для оптимизации параметров Солнечной системы. Но одно дело, что у них там написано, а другое дело то, что я вижу. Да и мой собственный опыт с коэффициентами в формуле Планка (1t) не очень показательный, т.к. оптимизировал я там только 3 коэффициента, а 4-ый фактор (температура излучения) мною принудительно задавался для повышения качества информации при проведение вычислительных экспериментов по почти D-оптимальному плану Бокса для четырех факторов. К тому же второй коэффициент (показатель степени при частоте излучения v) мне надо было не столько оптимизировать, сколько подтвердить, что он равен 3, как это следовало из формулы Вина. Таким образом, я по большому счету оптимизировал только 2 параметра (фактора) и, следовательно, у меня могли быть только двойные смешанные взаимодействия, а это прекрасно воспроизводится полиномом 2-ой степени (2), который я получаю после обработки данных вычислительных экспериментов.

А кто может ответить на вопрос есть ли в нашей системе, которую мы исследуем, смешанные взаимодействия выше парных, т.е. тройные или четверные. А может быть даже есть и не только линейные взаимодействия, но и квадратичные. К сожалению, ответить на эти вопросы никто не может. Да, наверное, никто не сможет ответить и на то, как это скажется на описание поверхности отклика при таких условиях. По этому, я на всякий случай (не очень доверяя всему, что написано в учебниках) решил провести маленькое исследование по оптимизации по критериям dY и dY^2 параметров простейших математических выражений, которые будут имитировать поведение различных систем. И первым делом я решил взять чуть ли не самый сложный случай с четверным взаимодействием, где вдобавок одно взаимодействие еще и не линейно, т.е. всем Вам известный закон тяготения Ньютона (3t) и попробовать оптимизировать его параметры. В принципе, мы можем с законом тяготения провести и натурные эксперименты. Правда, не с самим законом тяготения для масс, а с законом тяготения для зарядов (закон Кулона), где даже аналог гравитационной постоянной можем изменять, распологая различные диэлектрики между зарядами. Но речь сейчас идет не о том, можем ли мы воспроизвести эксперименты на реальном объекте или на его модели, а о том, можем ли мы, уже даже зная аналитическую формулу, отражающую отклик системы на наши воздействия на нее, чисто с математической точки зрения получить оптимальные значения системы по примененному мною критерию dY, т.е. по разнице между откликом системы и известным оптимальным значением.

Может возникнуть вопрос а зачем вообще надо проводить исследования для получения аппроксимации (2), если у нас уже есть аналитическая формула закона тяготения. А затем, что, мы сейчас просто проверяем на что способны методы многофакторного планирования, чтобы заранее знать, что от них ожидать. Ведь когда мы исследуем какую то сложную систему, то нам надо проводить натурные или вычислительные эксперименты, чтобы получить хотя бы уравнение регрессии (2), т.к. никакие аналитические выражения для критерия оптимизации при исследование самого объекта нам не известны вообще, а аналитическая формула, по которой вычисляется критерий оптимизации в моделях объекта, даже если и удастся такую получить в развернутом виде, может уместиться только на десятках или сотнях страниц, что делает ее не пригодной для аналитических методов оптимизации. А уравнение регрессии (2), т.е. полином 2-го порядка, который мы получаем при многофакторном планировании, очень удобен для этого и по этому мы и постараемся его получить по критерию dY для тестируемых систем. А т.к. в программе Solsys5 у меня по формуле (3) рассчитывается значение комплексного критерия оптимизации (целевой функции) в каждом из 24 экспериментов, а отклик системы в наших тестовых примерах определяется не по комплексному критерию, то мы можем, для оптимизации параметров по критерию dY в тестовых выражениях (3t10t), формулу (3) упростить до выражения (4)

Yu0(U) = Abs((YRas(U) - Yopt) / Yopt) (4)
Где: Yu0 (U) относительная разница между расчетным YRas(U) и оптимальным Yopt значением отклика системы, поведение которой имитирует одна из формул (3t10t), в U-ом эксперименте.

Результаты оптимизации параметров в формуле тяготения, по примененному мною критерию оптимизации dY, получились удручающие, т.к. аппроксимация критерия оптимизации, полученным уравнением регрессии (2), не лезла ни в какие ворота. Да Вы сами взгляните на полученные значения критерия оптимизации dY с использованием формулы закона тяготения (3t) и эти же значения по полученному уравнению регрессии (2) на нижеприведенном рисунке (верхняя часть рисунка), где маленькие синие кружки это критерий оптимизации рассчитанный с использованием формулы (3t) для выражения (3) в 24 экспериментах плана Бокса (номер соответствует абсциссе), а большие синие кружки это критерий оптимизации, рассчитанный по уравнению регрессии (2), для тех же значений параметров, что и в соответствующем эксперименте. При этом все факторы X1 X4, при выполнение плана Бокса, на нулевом уровне были равны единице, а интервалы их варьирования были 0,5, а значение Yopt бралось равным 1, т.е. я принимал, что оптимальные значения параметров X1 X4 в формуле (3t) равны 1 и получалось, что Yopt = 1*1*1/1=1. Для меня такой результат аппроксимации был большой неожиданностью, т.к. ранее в своей книге я сам закон тяготения, т.е. по критерию Y, аппроксимировал уравнением (2) и никаких проблем по качеству аппроксимации тогда не было. Я тут же аппроксимировал полученные значения Yu0(U) полиномом (2) и по критерию Y и выяснилось, что на этот раз у меня получились значительные погрешности в аппроксимации (синие маленькие и большие кружки на нижней части рисунка).

Yu = k0 + k1*X1 + k2*X2 + k3*X3 + k4*X4 +
+ k5*X1*X2 + k6*X1*X3 + k7*X1*X4 + k8*X2*X3 + k9*X2*X4 + k10*X3*X4 +
+ k11*X1^2 + k12*X2^2 + k13*X3^2 + k14*X4^2 (2)
где Yu (или dYu или dYu^2) критерий оптимизации, который надо минимизировать, X1 X4 - оптимизируемые параметры, а k0 k14 коэффициенты, которые мы получаем методом наименьших квадратов при статистической обработке значений Yu0(U) полученных в 24 экспериментах при разных значениях параметров X1 X4.

http://ser.t-k.ru/Ris/3t_Y_dY.gif (зеркало http://modsys.narod.ru/Ris/3t_Y_dY.gif)

Стал разбираться в чем дело и выяснил, что когда я ранее аппроксимировал закон притяжения, то задавал интервалы варьирования параметров примерно 20%, а сейчас задал их 50%. Уменьшил интервал варьирования до 20% и на сей раз опять получил приличный результат (черные кружки), т.е. получается, что все дело не в самом критерии оптимизации, а в интервалах варьирования, тем более, что при интервалах варьирования 80% (зеленые кружки), результат получился еще хуже, чем при 50%. И по критерию dY (в верхней части рисунка) тоже получается, что, чем больше интервал варьирования, тем хуже аппроксимация экспериментальных данных, но здесь уже не все так однозначно, т.к. получается, что и при интервале варьирования 20% аппроксимация экспериментальных данных получается не на много лучше. По этому я решил продолжить исследование с использованием других математических выражений имитирующих поведение системы.

YRas(U) = X1*v^X2* exp(-X3*v/X4) (1t)
YRas(U) = u (2t)
YRas(U) = X1 * X2 * X3 / X4^2 (3t)
YRas(U) = X1^2 * X2^2 (4t)
YRas(U) = X1 * X2 * X3 * X4 (5t)
YRas(U) = X1 + X1^2 (6t)
YRas(U) = X1 + X2 + X3 + X4 (7t)
YRas(U) = X1*X2 + X1*X3 + X1*X4 + X2*X3 + X2*X4 + X3*X4 (8t)
YRas(U) = X1^2 + X2^2 + X3^2 + X4^2 (9t)
YRas(U) = X1 + X2 + X3 + X4 + X1^2 + X2^2 + X3^2 + X4^2 +
+ X1*X2 + X1*X3 + X1*X4 + X2*X3 + X2*X4 + X3*X4 (10t)

При проведение вычислительных экспериментов по всем этим выражениям я принимал, что оптимальные значения всех параметров в этих формулах равны единице и находил сначала оптимальное значение отклика системы Yopt, а затем задавал в соответствие с планом различные значения параметров X1 X4 и, вычислив значение YRas(U), находил значение критерия оптимизации в U-ом эксперименте. Значение критерия Y я брал равным YRas(U), критерия dY я вычислял по формуле (4), а критерий dY^2 определял возведя критерий dY в квадрат. Затем, обработав данные по 24 значениям этих критериев Yu0(U), я получал коэффициенты k0 k14 для аппроксимации (2) по которой вычислял значение Yu для тех же значений параметров X1 X4, что были в 24 экспериментах по плану эксперимента, и сравнивал их с 24 значениями Yu0(U). А по результатам сравнения я выставлял оценку почти D-оптимальному плану Бокса по качеству аппроксимации экспериментальных данных выражением (2) по различным критериям оптимизации. Сравнение я проводил графически определяя попало ли значение Yu0(U), которое на приведенных выше рисунках отражено в 24 экспериментах маленькими кружками, внутрь большого кружка отражающего значение Yu. При этом графический масштаб для вывода данных выбирался так, чтобы все данные от минимального до максимального значения по ординате укладывались в интервале от 1 до 7 сантиметров. А оценки качеству аппроксимации значений критериев оптимизации Y, dY и dY^2, полученных по плану Бокса полиномом (2), я производил по пятибальной шкале, определяя наилучшую оценку по наихудшим результатам, а затем, полученные результаты, оформил в виде таблицы. Критерии оценок были такими

5 баллов все 24 маленьких кружка, т.е. значения Yu0(U), находятся внутри соответствующего большого кружка, центр которого соответствует ординате значения Yu. При этом диаметр большого кружка в два раза больше диаметра маленького.
4 балла - все 24 маленьких кружка или находятся внутри больших кружков или хотя бы касаются его с наружной стороны. При этом я также указываю в таблице, рядом с оценкой, в знаменателе количество экспериментов, по которым была выставлена эта наихудшая оценка. Так, если в 23 случаях маленькие кружки находятся внутри больших, а в одном случае маленький кружок пересекается с большим или касается его с наружной стороны, то будет указана оценка 4/1.
3 балла если все маленькие кружки находятся хотя бы на расстояние 1-го диаметра большого кружка от его окружности. При этом, если 20 маленьких кружков находятся внутри больших, 1 кружок пересекается с большой окружностью, а три маленьких кружка находятся на расстояние 1-го большого диаметра от его окружности, то оценка будет 3/3. Таким образом, в знаменателе указывается только количество кружков не прошедших по более высоким оценкам.
2 балла если все маленькие кружки находятся хотя бы на расстояние 2-х диаметров большого кружка от его окружности.
1 балл - если хотя бы один маленький кружок находится на расстояние более 2-х диаметров большого кружка от его окружности.

Таблица 15. Оценки качества аппроксимации полиномом 2-ой степени различных критериев оптимизации при проведение вычислительных экспериментов по почти D-оптимальному плану Бокса, когда эти критери расчитывались по разным математическим выражениям (имитирующим поведение системы) и при разных интервалах варьирования оптимизируемых факторов (параметров), которые на нулевых уровнях принимали значение 1.

Формула_________критерий_Y________критерий_dY________критерий_dY^2
Интервалы_+/-0,2_+/-0,5_+/-0,8_+/-0,2_+/-0,5_+/-0,8_+/-0,2_+/-0,5_+/- 0,8
3t__________4/2___2/1___1/2______3/3___2/2___1/2_____1/2___1/3___1/6
4t________ ___5____4/4___3/4______3/4___3/4___3/4_____2/4___2/4___2/4
5t___________5___3/10___1/1__ ____2/1___1/2___1/6_____1/2___1/6___1/6
6t___________5____5_____5________5_____5_____5__ ____5_____5_____5
7t___________5____5_____5______3/16___3/16___3/16____5_____5_____5
8t___________5____5_____5________2/1___2/1___2/1_____3/2___2/2___1/2
9t___________5___ _5_____5_______3/11___3/5___3/6_____4/8___4/9___3/6
10t__________5____5_____5______3/16_ __3/15___3/15____4/2___3/2___2/2

Что можно сказать о полученных данных. То, что в тестах 6t10t по критерию Y результаты будут отличными я и не сомневался, т.к. для аппроксимации именно таких взаимодействий полиномы 2-го порядка и приспособлены, но вот почему по критериям dY и dY^2 получились такие плачевные результаты это для меня не понятно. Ясно было и то, что в тестах 3t5t результаты будут хуже чем в тестах 6t10t, т.к. линейные взаимодействия выше парных в примененном нами плане Бокса смешаны с другими взаимодействиями и по этому их нельзя выделить (про парное квадратичное взаимодействие затрудняюсь сказать что то определенное), но и в этом случае я никак не ожидал, что по критериям dY и dY^2 будут такие плохие результаты. По этому остается только надеяться на то, что в исследуемой нами Солнечной системе тройных и выше взаимодействий не будет, а то нам придется очень долго мучиться пока мы найдем область оптимума. И я уже подумываю о включение в 6-ю версию программы ротатабельного плана, у которого информация более равномерно размазана по всем направлениям, т.е. одинаковая на всех направлениях на равных расстояниях от центра плана, а, учитывая то, что звездное плечо у этого плана равно 2, это позволит прощупать факторное пространство на ту же глубину при меньшем значении интервала варьирования для центральной части плана, и по этому этот план может оказаться нам очень полезен именно при поиске области оптимума.

Однако, давайте наконец то посмотрим насколько удачно наш полином (2) аппроксимирует смоделированные нами процессы в Солнечной системе по критериям dY и dY^2. На нижеприведенном рисунке (в верхней части) показаны значения Yu0(U) полученные после обработки данных вычислительных экспериментов на классической математической модели Солнечной системы с учетом скорости распространения гравитации по критерию dY по формуле (3) (маленькие синие кружки) и по критерию dY^2 (маленькие красные кружки) и соответствующие им значения Yu (большие кружки) полученные по уравнениям (2). Причем в уравнение (3) значения рассчитаны по комплексному критерию, т.е. с учетом вековых смещений всех J-х параметров для Меркурия, т.е. J=1 4, а I=1 и с весовыми коэффициентами kVesa(I, J)=100. Как видим по критерию dY аппроксимацию можно оценить на твердую тройку, а если бы не 23 и 24 эксперименты, где у нас был очень большой интервал варьирования скорости гравитации, то возможно бы получилась и четверка. Кстати и интервал варьирования скорости VZsys тоже великоват, т.к. в экспериментах 17 и 18, 19 и 20, 21 и 22, 23 и 24 разность значений критерия между этими парными экспериментами должна быть примерно равна. И по критерию dY^2 аппроксимацию тоже можно оценить на тройку, а если уменьшить интервалы варьирования VZsys и Vgr то тоже может получиться хорошая оценка.

http://ser.t-k.ru/Ris/Plan3.gif (зеркало http://modsys.narod.ru/Ris/Plan3.gif)

Таким образом, быстрее всего, в моделируемой нами системе нет взаимодействий выше парных и аппроксимацию в виде полинома 2-ой степени можно использовать для поиска оптимума. Кстати, если кому-то показалось, что методы многофакторного планирования не очень надежны, то хочу Вас разочаровать, т.к. в основе этих методов лежит очень мощный математический аппарат, хотя и зародились эти методы для исследований в сельском хозяйстве, т.е. в области не очень блещущей математическими изысками. А в доказательство хочу Вам показать маленький фокус, на котором я всегда тестирую все свои новые программы, где применяю многофакторное планирование. На нижней части второго рисунка Вы видите значения Yu0(U) по критерию Y для теста 2t и эти же значения полученные по уравнению регрессии (2) для Yu. А фокус заключается в том, что не зависимо от уровней варьирования факторов мы принимаем значение Yu0(U) в очередном эксперименте просто равным номеру этого эксперимента, т.е. получается, что уравнением (2) мы аппроксимировали почти полный хаос, но результат то получился вполне приличным.

Я конечно же понимаю, что некоторые болельщики за ту или иную теорию уже заждались, когда же я объявлю окончательный результат матча, а я тут рассказываю про скучные двойные взаимодействия и звездные плечи, а им гораздо интереснее было бы услышать про двойные звезды. Но, к сожалению, наука это не сказка и быстро только эти самые сказки сказываются, но не скоро дело делается, а т.к. вопрос мы рассматриваем очень сложный то до окончания этого матча пока еще очень далеко. Хотя, не скрою, что после выполнения третьего плана, я начал уже склоняться к тому, что классическая модель с учетом скорости распространения гравитации (если не в этом матче, то в этом эпизоде уж точно) проиграла. Правда, это, конечно же, не означает, что в этом эпизоде победила ОТОшная модель, т.к. я ее в этих условиях вообще еще не тестировал. А вот, учитывая проведенное мною исследование по описанию систем имитируемых формулами 3t10t, можно пока еще все трудности по нахождению оптимума объяснить именно трудностью описания поверхности отклика системы по критериям dY или dY^2. И хотя, конечно же, как многие догадались, я болею в этом матче именно за классическую модель, но я здесь в лучшем случае арбитр и по этому кто победит покажет только игра. А я сейчас сделаю описание к 5-ой версии программы Solsys и выложу ее для скачивания, чтобы самые нетерпеливые болельщики сами смогли провести на ней вычислительные эксперименты и попытаться найти оптимум для наших параметров системы (если, конечно же, он там есть для классической модели Солнечной системы).

С наилучшими пожеланиями Сергей Юдин.