Comments 33
А с табличными методами (интерполяция многочленом, например) сравнивали?
Ну или с более шустрыми методами численного поиска корня, типа метода Ньютона?
В таких случаях стоит применять более продвинутые инструменты.
Вот пример автоматической генерации текста функции на C для PCHIP сплайна на основе табулированных данных в среде MATLAB
Генератор исходников функции интерполяции PCHIP сплайном
max_arg = 180;
min_arg = -180;
file_name = "Speed_control_by_angle";
arr_postfix="speed_scale";
arg_name = "rod_angle";
spline_arr_name = sprintf('%s_%s', 'spline_arr',arr_postfix);
lookup_arr_name = sprintf('%s_%s', 'lookup_arr',arr_postfix);
max_lookup_arg_name = upper(sprintf('%s_%s', 'MAX_LOOKUP',arr_postfix));
min_lookup_arg_name = upper(sprintf('%s_%s', 'MIN_LOOKUP',arr_postfix));
arg_step_name = upper(sprintf('%s_%s', 'STEP_LOOKUP',arr_postfix));
% Табулированные входные данные
arg_step = 10;
x = [min_arg -170 -160 -150 -140 -130 -120 -110 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 max_arg];
y = [0.99340593 0.939399192 0.872323991 0.823953826 0.799539693 0.794322777 0.802430832 0.819113434 0.84086591 0.865140001 0.890064412 0.914242756 0.936618244 0.956384977 0.972929976 0.985795264 0.994653173 0.999290618 0.999599748 0.998966483 0.987305079 0.974992332 0.958945791 0.939603646 0.917554893 0.893575028 0.868680621 0.844212815 0.821964837 0.804373374 0.794787874 0.795421548 0.819088361 0.863978024 0.929558487 0.988495752 0.993411949];
pp = pchip(x,y);
xx = linspace(min_arg,max_arg,1000); % Выводим контрольный график на 1000 точек
plot(x,y,'o',xx,ppval(pp,xx),'-');
%pp.coefs
[m,n] = size(x); % n - количество элементов в массиве x
left_y = y(1);
right_y = y(n);
file_name_c = sprintf('%s.c', file_name);
file_name_h = sprintf('%s.h', file_name);
fid = fopen(file_name_c,'w','n','windows-1251');
fprintf(fid,'#include "App.h"\n');
fprintf(fid,'\n\n');
fprintf(fid,'typedef struct \n');
fprintf(fid,'{\n');
fprintf(fid,' float arg;\n');
fprintf(fid,' float k3;\n');
fprintf(fid,' float k2;\n');
fprintf(fid,' float k1;\n');
fprintf(fid,' float k0;\n');
fprintf(fid,'} T_splane_rec;\n');
fprintf(fid,'\n');
fprintf(fid,'\n');
fprintf(fid,'typedef struct \n');
fprintf(fid,'{\n');
fprintf(fid,' uint32_t indx;\n');
fprintf(fid,'} T_sline_lookup_rec;\n');
fprintf(fid,'\n');
fprintf(fid,'\n');
fprintf(fid,'#define %s %.1ff\n',max_lookup_arg_name, max_arg);
fprintf(fid,'#define %s %.1ff\n',min_lookup_arg_name, min_arg);
fprintf(fid,'#define %s %.0fUL\n' ,arg_step_name, arg_step);
fprintf(fid,'\n');
fprintf(fid,'// Spline PCHIP for %s\n', arr_postfix);
fprintf(fid,'// Input x vector = [%s]\n', join(pad(string(x),6),', '));
fprintf(fid,'// Input y vector = [%s]\n', join(pad(string(y),6),', '));
% Печатаем массив коэффициентов сплайна
fprintf(fid,'\n');
fprintf(fid,'T_splane_rec %s[%d] = \n',spline_arr_name, n-1);
fprintf(fid,'{\n');
for i = 1:(n-1)
s = sprintf(' {%4d, %14.5ef, %14.5ef, %14.5ef, %14.5ef}, \n',x(i), pp.coefs(i,1), pp.coefs(i,2), pp.coefs(i,3), pp.coefs(i,4));
fwrite(fid, s);
end
fprintf(fid,'};\n');
% Печатаем массив соответствия угла индексу в предыдущей таблице
angle = min_arg:arg_step:max_arg;
[l,k] = size(angle); % k - количество элементов в массиве angle
fprintf(fid,'\n');
fprintf(fid,'T_sline_lookup_rec %s[%d] = \n',lookup_arr_name, k-1);
fprintf(fid,'{\n');
indx = 1;
for i = 1:k
if angle(i)>=x(indx+1)
indx = indx + 1;
if indx > (n -1)
break;
end
end
s = sprintf(' {%4d }, // %4d \n', indx-1, angle(i));
fwrite(fid, s);
end
fprintf(fid,'};\n');
fprintf(fid,'/*-----------------------------------------------------------------------------------------------------\n');
fprintf(fid,' Нелинейная функция на основе кубических PCHIP сплайнов \n');
fprintf(fid,' https://www.mathworks.com/help/matlab/ref/pchip.html \n');
fprintf(fid,'-----------------------------------------------------------------------------------------------------*/\n');
fprintf(fid,'float Calk_spline_%s( float %s)\n',arr_postfix, arg_name );
fprintf(fid,'{\n');
fprintf(fid,' uint32_t indx;\n');
fprintf(fid,' float f;\n');
fprintf(fid,'\n');
fprintf(fid,' if (%s >= %s)\n', arg_name, max_lookup_arg_name );
fprintf(fid,' {\n');
fprintf(fid,' return %ef;\n', right_y);
fprintf(fid,' }\n');
fprintf(fid,' if (%s < %s)\n', arg_name, min_lookup_arg_name);
fprintf(fid,' {\n');
fprintf(fid,' return %ef;\n', left_y);
fprintf(fid,' }\n');
fprintf(fid,' indx = (uint32_t)(%s - %s);\n', arg_name, min_lookup_arg_name);
fprintf(fid,' indx = %s[indx/%s].indx;\n', lookup_arr_name, arg_step_name);
fprintf(fid,' %s = %s - %s[indx].arg;\n', arg_name, arg_name, spline_arr_name);
fprintf(fid,'\n');
fprintf(fid,' f = %s ', arg_name);
fprintf(fid,'* ((%s[indx].k3 ', spline_arr_name);
fprintf(fid,'* %s ', arg_name);
fprintf(fid,'+ %s[indx].k2) ', spline_arr_name);
fprintf(fid,'* %s ',arg_name);
fprintf(fid,'+ %s[indx].k1) ', spline_arr_name);
fprintf(fid,'+ %s[indx].k0;\n', spline_arr_name);
fprintf(fid,'\n');
fprintf(fid,' return f;\n');
fprintf(fid,'}\n');
fclose(fid);
fid = fopen(file_name_h,'w','n','windows-1251');
fprintf(fid,'#ifndef %s\n',upper(file_name));
fprintf(fid,'#define %s\n\n',upper(file_name));
fprintf(fid,'float Calk_spline_%s( float %s);\n\n',arr_postfix, arg_name );
fprintf(fid,'#endif // %s\n',upper(file_name));
fclose(fid);И на выходе получаем такую функцию
Получившийся исходник функции интерполяции на языке C
#include "App.h"
typedef struct
{
float arg;
float k3;
float k2;
float k1;
float k0;
} T_splane_rec;
typedef struct
{
uint32_t indx;
} T_sline_lookup_rec;
#define MAX_LOOKUP_DOOR_SPEED_SCALE 180.0f
#define MIN_LOOKUP_DOOR_SPEED_SCALE -180.0f
#define STEP_LOOKUP_DOOR_SPEED_SCALE 10UL
// Spline PCHIP for door_speed_scale
// Input x vector = [-180 , -170 , -160 , -150 , -140 , -130 , -120 , -110 , -100 , -90 , -80 , -70 , -60 , -50 , -40 , -30 , -20 , -10 , 0 , 10 , 20 , 30 , 40 , 50 , 60 , 70 , 80 , 90 , 100 , 110 , 120 , 130 , 140 , 150 , 160 , 170 , 180 ]
// Input y vector = [0.99341, 0.9394, 0.87232, 0.82395, 0.79954, 0.79432, 0.80243, 0.81911, 0.84087, 0.86514, 0.89006, 0.91424, 0.93662, 0.95638, 0.97293, 0.9858, 0.99465, 0.99929, 0.9996, 0.99897, 0.98731, 0.97499, 0.95895, 0.9396, 0.91755, 0.89358, 0.86868, 0.84421, 0.82196, 0.80437, 0.79479, 0.79542, 0.81909, 0.86398, 0.92956, 0.9885, 0.99341]
T_splane_rec spline_arr_door_speed_scale[36] =
{
{-180, 7.05244e-07f, -7.23948e-05f, -4.74725e-03f, 9.93406e-01f},
{-170, 1.81073e-05f, -2.53468e-04f, -5.98357e-03f, 9.39399e-01f},
{-160, 8.08325e-06f, -2.46076e-06f, -5.62073e-03f, 8.72324e-01f},
{-150, 7.78170e-06f, 2.53899e-06f, -3.24497e-03f, 8.23954e-01f},
{-140, 1.83701e-06f, 1.54290e-05f, -8.59683e-04f, 7.99540e-01f},
{-130, -5.30366e-06f, 1.34117e-04f, 0.00000e+00f, 7.94323e-01f},
{-120, -3.56959e-06f, 9.33974e-05f, 1.09125e-03f, 8.02431e-01f},
{-110, -1.67758e-06f, 4.54690e-05f, 1.88832e-03f, 8.19113e-01f},
{-100, -1.00902e-06f, 2.33890e-05f, 2.29442e-03f, 8.40866e-01f},
{ -90, -7.08159e-07f, 1.03762e-05f, 2.45950e-03f, 8.65140e-01f},
{ -80, -5.68971e-07f, 2.01606e-06f, 2.45457e-03f, 8.90064e-01f},
{ -70, -5.18604e-07f, -3.47915e-06f, 2.32420e-03f, 9.14243e-01f},
{ -60, -5.30158e-07f, -6.93474e-06f, 2.09904e-03f, 9.36618e-01f},
{ -50, -6.02108e-07f, -8.65836e-06f, 1.80129e-03f, 9.56385e-01f},
{ -40, -7.63660e-07f, -8.45999e-06f, 1.44749e-03f, 9.72930e-01f},
{ -30, -1.13611e-06f, -4.97945e-06f, 1.04920e-03f, 9.85795e-01f},
{ -20, -2.60753e-06f, 1.15724e-05f, 6.08773e-04f, 9.94653e-01f},
{ -10, -3.86374e-08f, -2.31855e-06f, 5.79623e-05f, 9.99291e-01f},
{ 0, 6.52355e-08f, -6.98501e-06f, 0.00000e+00f, 9.99600e-01f},
{ 10, 1.01433e-05f, -2.06034e-04f, -1.20129e-04f, 9.98966e-01f},
{ 20, -1.28658e-06f, 9.52062e-06f, -1.19782e-03f, 9.87305e-01f},
{ 30, 6.18344e-07f, -2.73104e-05f, -1.39338e-03f, 9.74992e-01f},
{ 40, 5.36445e-07f, -2.33770e-05f, -1.75409e-03f, 9.58946e-01f},
{ 50, 5.16752e-07f, -1.95855e-05f, -2.06070e-03f, 9.39604e-01f},
{ 60, 5.57351e-07f, -1.56342e-05f, -2.29738e-03f, 9.17555e-01f},
{ 70, 6.80971e-07f, -1.14680e-05f, -2.44286e-03f, 8.93575e-01f},
{ 80, 9.51197e-07f, -7.39740e-06f, -2.46793e-03f, 8.68681e-01f},
{ 90, 1.54322e-06f, -4.86042e-06f, -2.33052e-03f, 8.44213e-01f},
{ 100, 3.12608e-06f, -1.06995e-05f, -1.96476e-03f, 8.21965e-01f},
{ 110, 6.76174e-06f, -3.93798e-05f, -1.24093e-03f, 8.04373e-01f},
{ 120, -3.30481e-08f, 6.66722e-06f, 0.00000e+00f, 7.94788e-01f},
{ 130, -1.51060e-05f, 3.75386e-04f, 1.23430e-04f, 7.95422e-01f},
{ 140, -5.48865e-06f, 1.93850e-04f, 3.09933e-03f, 8.19088e-01f},
{ 150, -1.57819e-05f, 2.80650e-04f, 5.32974e-03f, 8.63978e-01f},
{ 160, -4.67175e-05f, 4.35731e-04f, 6.20817e-03f, 9.29558e-01f},
{ 170, -7.57014e-07f, -3.40217e-05f, 9.07538e-04f, 9.88496e-01f},
};
T_sline_lookup_rec lookup_arr_door_speed_scale[36] =
{
{ 0 }, // -180
{ 1 }, // -170
{ 2 }, // -160
{ 3 }, // -150
{ 4 }, // -140
{ 5 }, // -130
{ 6 }, // -120
{ 7 }, // -110
{ 8 }, // -100
{ 9 }, // -90
{ 10 }, // -80
{ 11 }, // -70
{ 12 }, // -60
{ 13 }, // -50
{ 14 }, // -40
{ 15 }, // -30
{ 16 }, // -20
{ 17 }, // -10
{ 18 }, // 0
{ 19 }, // 10
{ 20 }, // 20
{ 21 }, // 30
{ 22 }, // 40
{ 23 }, // 50
{ 24 }, // 60
{ 25 }, // 70
{ 26 }, // 80
{ 27 }, // 90
{ 28 }, // 100
{ 29 }, // 110
{ 30 }, // 120
{ 31 }, // 130
{ 32 }, // 140
{ 33 }, // 150
{ 34 }, // 160
{ 35 }, // 170
};
/*-----------------------------------------------------------------------------------------------------
Нелинейная функция на основе кубических PCHIP сплайнов
https://www.mathworks.com/help/matlab/ref/pchip.html
-----------------------------------------------------------------------------------------------------*/
float Calk_spline_door_speed_scale( float rod_angle)
{
uint32_t indx;
float f;
if (rod_angle >= MAX_LOOKUP_DOOR_SPEED_SCALE)
{
return 9.934119e-01f;
}
if (rod_angle < MIN_LOOKUP_DOOR_SPEED_SCALE)
{
return 9.934059e-01f;
}
indx = (uint32_t)(rod_angle - MIN_LOOKUP_DOOR_SPEED_SCALE);
indx = lookup_arr_door_speed_scale[indx/STEP_LOOKUP_DOOR_SPEED_SCALE].indx;
rod_angle = rod_angle - spline_arr_door_speed_scale[indx].arg;
f = rod_angle * ((spline_arr_door_speed_scale[indx].k3 * rod_angle + spline_arr_door_speed_scale[indx].k2) * rod_angle + spline_arr_door_speed_scale[indx].k1) + spline_arr_door_speed_scale[indx].k0;
return f;
} PCHIP сплайны хороши тем, что не уходят на краях неизвестно куда, а уходят по прямой касательной к сплайну на конечной точке.
Генератор, а не шаблон здесь выбран для того чтобы обеспечить глобальную уникальность имён внутри функции, скорость и прозрачность для отладки.
Если я правильно понял, то исходное уравнение (1) при заданных значениях delta и Day для неизвестной phi аналитически сводится к квадратному уравнению относительно неизвестного Sin(phi) .Для этого заменяем синус через косинус и обозначаем синус как Х.
В результате получается квадратное уравнение.
Решение в виде формулы известно, от которого арксинус.
На Си это будет как два пальца и без численных методов.
Попробуйте вывести конечную формулу. Я запутался когда пробовал.
Это же функция от двух переменных, и обратная функция тоже будет от двух переменных. Аналитическое решение вроде бы существует, но у меня в Вольфраме оно получилось чересчур громоздким - вычислительно точно будет ничуть не легче. Аппроксимировать функцию многих переменных тоже удовольствие такое себе. Поэтому автор сделал вполне разумный выбор.
Возможно не понял задачу, но я выше написал , как я ее понял. повторю:
Необходимо найти phi из уравнения (1) при заданных значениях delta и Day. Или не так?
Ну, delta и Day - это два аргумента, а не один. Ниже я привёл своё аналитическое решение, с численным от автора он совпадает. Если у вас есть другое, более простое решение - будет очень интересно посмотреть.
Необходимо найти phi из уравнения (1) при заданных значениях delta и Day. Или не так?
Да. Именно так.
Тогда это уравнение с одним неизвестным. Заменяем Sin(phi) через Cos(phi) а Cos(phi) обозначаем X. Получаем квадратное уравнение ( X в степени 2) Что сложного?
------------------
A=COS(PI*DAY/24)
B=-SIN(51')
C=A*Cos(delta)
E=Sin(delta)
C*COS(phi)=B-E*SIN(phi)
COS(phi)^2=1-Sin(phi)^2
X=Sin(phi)
C^2*(1-SIN(phi)^2)=B^2- 2*B*E*Sin(phi) +Sin(phi)^2
C^2*(1-X^2)=B^2-2*B*E*X+X^2
Я понял так, что солнечное склонение тоже известно. Автор пишет что оно меняется раз в день.
В итоге не одна, а две константы и уравнение с одним неизвестным phi. Что не так?
Что не так? Что это константа присутствует и в числителе, и в знаменателе, и с аргументом связана нелинейно - поэтому её нельзя просто так взять и выкинуть, чтобы упростить уравнение для его решения. Квадратное уравнение здесь никак не получается, это даже по графику видно.
Давайте уточним, а то что-то не сходится.
Мы говорим с Вами о формуле (1), верно? Эта формула от одной переменной (см слева в скобках) В этой формуле задано delta и Day .
Если мы все константы, в том числе и синусы и косинусы от них обозначим буквами А В С,
то в формуле слева останутся лишь синус и косинус от переменной phi.
косинус заменяем на корень квадратный из 1 минус синус в квадрате.
Далее син phi обозначаем Х. В результате получаем квадратное уравнение от X.
Что не так?
Вы пишете что расчёт по формуле с малым шагом 0.001 градус будет долго выполняться на микроконтроллере. Но ведь можно рассчитать таблицу заранее и обращение к ней по сути - одна машинная команда. Или в используемом микропроцессоре мало ПЗУ ?
Визуально график достаточно простой, возможно, его можно интерполировать сплайнами.
Про точность, приведённую в статье в 0.001 градус достаточно float. Тип double только 'жирные" микроконтроллеры поддерживают аппаратно.
Вообще-то для подобных задач полиномы Чебышева придуманы.
Аналитическое решение:
Может, вручную можно ещё подсократить.
phi[day,delta] = -ArcTan[(0.7771459614569708 + Cos[0.1308996938995747*day])*Cot[delta]];
А чем вас Ньютоновской метод решения нелинейных уравнений не устроил? У вас случай идеальный, функция монотонная, начальное приближение есть, производную можно численно посчитать совершенно примитивными методами, за 3-4 итерации можно добиться гарантированной сходимости с очень высокой точностью в вашем случае.
26-го декабря у Солнца склонение не -18.74, а -23 с чем-то градуса. Почти как в солнцестояние.
Спасибо за статью. Задача ваша интересна с научной точки зрения, а слишком прямолинеен. Возникает вопрос - зачем проводить так много вычислений на конечном устройстве? Можно ж подумать заранее и все вычислить при разработке, а для МК останется расчет на порядки проще.
Формализуем задачу: у вас есть Day(phi, delta) заданная вот той формулой с тригонометрией и ее надо преобразовать к phi(Day, delta). Само преобразование - решение трансцендентного уравнения которое вычислительно сложно, нужно вычислять функцию с тригонометрией да еще и несколько раз. Я бы предложил аппроксимировать решение дробно рациональным полиномом или еще более сложной функцией. Аппроксимацию можно вести методом наименьших квадратов по большому массиву заранее расчитанных точек. А вот по подбору аппроксимирующих функций у нас полная свобода выбора.
При этом множество точек на которых мы проводим аппрксимацию нужно строить с учетом наличия особенностей - полярного дня и ночи. Я предлагаю такой вариант: решить вспомогательную задачу нахождения Phi24(delta)=phi(24, delta) и Phi0(delta)=phi(0,delta), т.е. широт полярного дня и ночи
Итоговое решение можно представить в виде phi(Day, delta)=Phi24(delta)*F(Day, delta) + Phi0(delta) (1-F(Day, delta)) где F некий рациональный полином равный 1 при Day=24 и 0 при Day=0. Вид F еще стоит конструировать, но тут без экспериментов не скажу. Такое представление должно давать хорошие результаты вблизи зон дня и ночи благодаря особенности построения, функций а более менее правдивый результат при низких широтах даст подбор коэффициентов. В итоге вместо тригонометрии и кучи итераци будет вычисление трех рациональных полиномов. Это в разы быстрее.
Если не интересует точность вблизи зон солнцестояния, можно вообще подобрать один рациональный полиним будет еще быстрее вычисляться. Заморочка с факторизацией описанная выше просто обусловлена фактом что внутренняя структура аппрксимаций должна учитывать поведение вблизи особых точек. А в середине интервала где все гладко и без особенностей любая аппроксимация справится.
На мой вкус вычисления стоит проводить в обезразмеренных величинах, т.е. delta обезразмерить до [-1..1] вместо [-23...23] градусов, а Day отобразить на [0..1] поделив на 24, но это вкусовщина
Любопытно было бы получить обратный расчет необходимой точности прибилжения исходя из качества входных данных. Вполне возможно что для данного приложения в целевой точности в 0.001° окажутся лишние нули. В таком случае и спектр применимых методов может расширится.
Любопытно было бы получить обратный расчет необходимой точности приближения исходя из качества входных данных.
Это уже проделано проделано. Вот в этом тексте https://habr.com/ru/articles/687640/
Формулы (1.2) (1.4)
Зачем Программисту Микроконтроллеров Численные Методы?