Как построить задачу поиска оптимальной точки?

Question

[ivandzemianchyk]

Есть 3 точки с коордниатами долготы и широты. Есть радиус окружности вокруг этих точек. Нужно построить модель для поиска точки, которая найдется в месте покрытия трех окружностей.

%[longtitude, Latitude ]
A = [17 24 18; 53 50 58];
B = [17 27 21; 53 56 45];
C = [17 43 24; 53 52 40];
radius = [6.53; 9.51; 14.76];

пробовал решить таким способом:

[x_opt fval exitflag output] = fminunc(@optFunction, [30; 30])

где optFunction.m, это

function f = optFunction(x)

%[dlugosc/ longtitude, szerokosc/Latitude ]
A = [17 24 18; 53 50 58];
B = [17 27 21; 53 56 45];
C = [17 43 24; 53 52 40];
d = [6.53; 9.51; 14.76];

calcVec = [1 1/60 1/3600];

latlon1 = [ A(2,:)*calcVec' A(1,:)*calcVec'];
latlon2 = [ B(2,:)*calcVec' B(1,:)*calcVec'];
latlon3 = [ C(2,:)*calcVec' C(1,:)*calcVec'];

f = (r1(latlon1, x, d))^2 + (r2(latlon2, x, d))^2 + (r3(latlon3, x, d))^2;
end

    function r1x = r1(latlon1, latlonx, d)
        r1x = (lldistkm(latlon1, latlonx))^2 - (d(1))^2;
    end

    function r2x = r2(latlon2, latlonx, d)
        r2x = (lldistkm(latlon2, latlonx))^2 - (d(2))^2;
    end
    
    function r3x = r3(latlon3, latlonx, d)
        r3x = (lldistkm(latlon3, latlonx))^2 - (d(3))^2;
    end

а файл lldistkm.m, это скачаный файл с www.mathworks.com/matlabcentral/fileexchange/38812...
а его код:

function [d1km d2km]=lldistkm(latlon1,latlon2)
% format: [d1km d2km]=lldistkm(latlon1,latlon2)
% Distance:
% d1km: distance in km based on Haversine formula
% (Haversine: http://en.wikipedia.org/wiki/Haversine_formula)
% d2km: distance in km based on Pythagoras’ theorem
% (see: http://en.wikipedia.org/wiki/Pythagorean_theorem)
% After:
% http://www.movable-type.co.uk/scripts/latlong.html
%
% --Inputs:
%   latlon1: latlon of origin point [lat lon]
%   latlon2: latlon of destination point [lat lon]
%
% --Outputs:
%   d1km: distance calculated by Haversine formula
%   d2km: distance calculated based on Pythagoran theorem
%
% --Example 1, short distance:
%   latlon1=[-43 172];
%   latlon2=[-44  171];
%   [d1km d2km]=distance(latlon1,latlon2)
%   d1km =
%           137.365669065197 (km)
%   d2km =
%           137.368179013869 (km)
%   %d1km approximately equal to d2km
%
% --Example 2, longer distance:
%   latlon1=[-43 172];
%   latlon2=[20  -108];
%   [d1km d2km]=distance(latlon1,latlon2)
%   d1km =
%           10734.8931427602 (km)
%   d2km =
%           31303.4535270825 (km)
%   d1km is significantly different from d2km (d2km is not able to work
%   for longer distances).
%
% First version: 15 Jan 2012
% Updated: 17 June 2012
%--------------------------------------------------------------------------

radius=6371;
lat1=latlon1(1)*pi/180;
lat2=latlon2(1)*pi/180;
lon1=latlon1(2)*pi/180;
lon2=latlon2(2)*pi/180;
deltaLat=lat2-lat1;
deltaLon=lon2-lon1;
a=sin((deltaLat)/2)^2 + cos(lat1)*cos(lat2) * sin(deltaLon/2)^2;
c=2*atan2(sqrt(a),sqrt(1-a));
d1km=radius*c;    %Haversine distance

x=deltaLon*cos((lat1+lat2)/2);
y=deltaLat;
d2km=radius*sqrt(x*x + y*y); %Pythagoran distance

end

Answer 1

[Mrrl]

Я бы её решал как точку пересечения трёх сфер в пространстве с последующей проекцией на плоскость ABC, а потом на поверхность Земли:
Пусть радиусы равны RA, RB, RC, а расстояния между точками - AB,AC,BC.
Вычислим величины
u=((RA^2-RB^2)/AB^2+1)/2
v=((RA^2-RC^2)/AC^2+1)/2
Это будут проекции точки пересечения сфер на отрезки AB и AC, выраженные в длинах этих отрезков.
Пусть a - угол в вершине А треугольника АВС (найдём по теореме косинусов).
Если Q - проекция точки пересечения сфер на плоскость ABC, то
AQ=(AB*(u-v*cos(a))+AC*(v-u*cos(a)))/(sin(a)^2),
т.е. Q=A*(1-(u+v)/(1+cos(a)))+B*((u-v*cos(a))/sin(a)^2)+C*((v-u*cos(a))/sin(a)^2),
где A,B,C - трёхмерные координаты точек. Тогда Q*(R/|OQ|) - искомая точка (|OQ| - расстояние от Q до центра Земли, R - радиус Земли).
Формулы не проверялись :(