План выполнил, данные обработал, а сейчас крепко задумался. К сожалению никакого глобального оптимума достичь не получилось, хотя данные для подобных планов многофакторного планирования получились просто отличные, если не считать двух последних экспериментов. Да Вы взгляните сами на результаты вычислительных экспериментов для Меркурия по плану Бокса, где факторы варьировались в следующих интервалах VXsys= 0+/-100, VYsys= 0+/-100, VZsys= 0+/-100, Vgr= 50*Vsv+/-30*Vsv (скорости системы по осям в км/с, а скорость гравитации в долях скорости света). Обозначения для вековых смещений параметров орбит те же, что были и раньше dAlfaP перигелий, dAlfaU узел восхождения, dBetta наклон орбиты и dEks эксцентриситет. Полученные смещения углов у меня даны в угловых секндах, а безразмерные значения вековых смещений эксцентриситетов увеличены мною в 3600 раз. В нулевой строке даны вековые смещения параметров орбит определенные по данным DE405 при неизвестных нам значениях факторов, которые и надо найти как оптимальные.

________dAlfaP__dAlfaU__dBetta__dEks
0__????___572,2__-450,0__- 21,4__0,074
1__++++__451,3__-1402__-230,5__0,208
2__+++-__234,8__-449,1__- 838,2__0,800
3__++-+__472,7__428,9___183,6__0,257
4__++--__ 273,5___2765___794,0__0,654
5__+-++__566,7__-1560__-263,2__0,302
6__+-+-__730,5___- 5206__-956,1__0,833
7__+--+__587,9___283,7__149,1___0,350
8__+--- ___791,0__2254___647,1___0,977
9__-+++__470,2__-1228__-194,1__-0,205
10_-++- ___289,0__-3858__-723,6__-1,071
11_-+-+___491,5__609,9___225,8__-0,154
12_-+-- ___376,5___3518___1002__-0,850
13_--++__586,1__-1386___-227,5__-0,110
14_--+- ___760,1__-4581__-848,4__-0,667
15_--+-___607,3__464,0___190,6__-0,061
16_---- ____845,2__3002___845,7__-0,468
17_+000__513,2__-591,8__-52,3____0,401
18_- 000__544,7__-307,6___10,0___-0,258
19_0+00__436,8__-329,4___5,8____-0,001
20_0- 00__621,4__-572,2___-48,6___0,148
21_00+0__512,3__-1984__-352,8___0,036
22_00- 0__545,9___963,4__312,8___0,114
23_000+__529,0__-450,9__-21,5___0,074
24_000- __529,1__-450,5___-21,5___0,074

Графики влияния каждого из факторов на критерий оптимизации по уравнению (2), когда остальные факторы зафиксированы на нулевых уровнях, я не привожу, т.к. по всем факторам они будут примерно такие же как и в первом плане для которого я их приводил. И определенный вывод по ним можно сделать только один глобальный оптимум по скорости гравитации находится где то за верхним уровнем, т.е. больше 80*Vsv. И я уже даже предчувствую критику знатоков планирования экспериментов за то, что я сразу стал ставить полные планы второго порядка, вместо того чтобы совершить крутое восхождение в область оптимума и уже там ставить полный план второго порядка, т.е. так как надо по классическим правилам планирования экспериментов. На что я могу ответить, что, во-первых, как я уже писал, мне хотелось побыстрее опробовать код программы по многофакторному планированию, а не излагать классическую теорию многофакторного планирования, во-вторых, некоторые поисковые эксперименты я все же проводил, и в третьих, самое главное, проводя вычислительные эксперименты на моделях систем, я всегда, вопреки рекомендациям классиков, использую только планы второго порядка.

Объясняется это тем, что, во-первых, число экспериментов, если мы не используем дробные реплики со смешанными коэффициентами (а дробные реплики отрицательно влияют на качество получаемой информации) при построение матрицы линейного плана для осуществления крутого восхождения, получается не на много меньше, чем при планах второго порядка (16 вместо 24), что при проведение вычислительных экспериментов не на много увеличит время их проведения (материальных затрат при этом практически никаких), а, во- вторых, чтобы методом крутого восхождения достичь области близкой к оптимуму надо выполнить 510 планов первого порядка и при этом общее число экспериментов может быть даже больше, т.к. информации от уравнения регрессии второго порядка можно получить гораздо больше и, следовательно, за меньшее число шагов достичь области оптимума. А кроме этого, если мы используем планы второго порядка, как при поисковых экспериментах, так и при изучении поверхности отклика в области оптимума, обработка данных производится по одной методике, что значительно упрощает код программы. Вот если бы эксперименты у нас были не вычислительные, т.е. не на математической модели, а натурные и нам в каждом эксперименте приходилось бы, например, запускать натуральную ракету к одной из планет, то тут можно было бы и подумать о различных градиентных методах для поиска области оптимума, т.к. после этих 48 экспериментов все исследования, наверное, и закончились бы, т.к. просто кончились бы деньги, выделенные на эти исследования. В нашем же случае, т.к. персональный компьютер потребляет не много электроэнергии, приходится жалеть только о затраченном времени.

А вообще-то, честно говоря, после этой серии экспериментов, я нахожусь в некотором замешательстве (неопределенности). Нечто подобное у меня было, когда я оптимизировал коэффициенты в формуле Планка, когда мне вследствие неопределенности решений никак не удавалось найти оптимум, но тогда я все таки нашел решение этого вопроса (кому интересно могу посмотреть это здесь http://ser.t-k.ru/Stat/Plank/plank.html зеркало http://modsys.narod.ru/Stat/Plank/plank.h tml ). А вот в этом случае я имею не только некоторую неопределенность, т.к., например, как видно из полученных данных, оптимум по всем вековым смещениям, кроме перигелия, может быть найден даже в уже заданных интервалах варьирования (кроме скорости гравитации). Вот только может получиться так, что каждый раз, выходя из разных точек гиперпространства, я буду каждый раз находить разные значения факторов для глобального оптимума. Но особенно меня смущает то, что пока выходит, что оптимальная скорость распространения гравитации лежит даже за пределами восьмидесяти скоростей света. Да, из литературных источников я знаю, что у многих исследователей теоретические значения скорости распространения гравитации получались очень большие, но одно дело об этом просто читать, а другое дело столкнутся с этим практически. Так что здесь у меня вдобавок к неопределенности существует еще и психологический барьер. Ведь мы уже все привыкли к тому, что скорость распространения для любых явлений Природы не может превышать скорость света.

Так что это действительно редкий случай, когда я принципиально не знаю что делать. А пока я буду осмысливать сложившуюся ситуацию, я решил выполнить еще один план, чтобы компьютер не простаивал, где решил задать на нулевом уровне скорость распространения гравитации 150 скоростей света. На всякий случай, решил немного изменить нулевые уровни и скоростей системы, а то не очень хорошо получается, когда на нулевых уровнях они все равны нулю. В общем, решил выполнить план с факторами VXsys= -50+/-100, VYsys= 0+/-100, VZsys= -50+/-100, Vgr= 150*Vsv+/-50*Vsv. Правда, в конечном итоге может оказаться, что все мои проблемы связаны с тем, что модель Солнечной системы, построенная на законах Ньютона (за исключением третьего закона), даже с учетом скорости распространения гравитации, просто не адекватно описывает процессы, происходящие в Природе, или, как говорят классики многофакторного планирования, у нас имеется не учтенный фактор, но пока у меня нет достаточных оснований, чтобы так заявить. А не учтенным фактором может быть и искривление пространства-времени, которое учитывается в ОТО, и динамическое давление гравитационного поля, которое учитывал Гербер, или еще что ни будь. Так что, как говорил Ю.Никулин //Будем искать//.

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

Обработал данные третьего плана многофакторного планирования и получил опять результаты, которые позволяют сделать только один вывод скорость распространения гравитации должна быть больше чем 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 кружок пересекается с большой окружностью, а три маленьких кружка находятся на расстояние