某货运公司为其所有卡车都配备了接收全球定位系统GPS信息的接收器，以此来确定卡车的位置。该系统约有30颗绕地球运行的卫星，每隔几秒钟，每颗卫星都同步地发出表明其准确位置及时间的信号。卡车上的接收器接听其中4颗卫星的信号，根据信号发出和到达的时间就能计算出卡车到卫星的距离（时间乘信号传播速度即光速）。假设4颗卫星的位置分别为（ai, bi, ci）(i=1,2,3,4)，卡车到这4颗卫星的距离分别为d1,d2,d3,d4.如何确定卡车的位置(x, y, z). 建立数学模型并给出求解模型的MATLAB程序。 ‎ ‌ ‎ ‌可以针对不同具体数据对模型或算法进行测试。 ‎ ‌比如取4颗同步卫星的位置为： ‎ ‌（15600,7540,20140）,（918760,2750,18610），（17610,14630,13480），（19170,610，18390）单位为千米， ‎ ‌卡车到这4颗卫星的距离分别为：22167.124千米，22604.821千米,24013.845千米,22670.775千米，求出该卡车的位置。 ‎

假设卡车的位置为(x, y, z)，则可以列出以下4个方程：

(x-a1)^2 + (y-b1)^2 + (z-c1)^2 = d1^2
(x-a2)^2 + (y-b2)^2 + (z-c2)^2 = d2^2
(x-a3)^2 + (y-b3)^2 + (z-c3)^2 = d3^2
(x-a4)^2 + (y-b4)^2 + (z-c4)^2 = d4^2

将每个方程展开并整理得到：

2xa1 - 2xa2 = d2^2 - d1^2 - a2^2 + a1^2 - b2^2 + b1^2 - c2^2 + c1^2
2ya1 - 2ya2 = d2^2 - d1^2 - a2^2 + a1^2 + b1^2 - b2^2 - c2^2 + c1^2
2za1 - 2za2 = d2^2 - d1^2 - a2^2 + a1^2 - b2^2 + b1^2 + c1^2 - c2^2

2xa1 - 2xa3 = d3^2 - d1^2 - a3^2 + a1^2 - b3^2 + b1^2 - c3^2 + c1^2
2ya1 - 2ya3 = d3^2 - d1^2 - a3^2 + a1^2 + b1^2 - b3^2 - c3^2 + c1^2
2za1 - 2za3 = d3^2 - d1^2 - a3^2 + a1^2 - b3^2 + b1^2 + c1^2 - c3^2

2xa1 - 2xa4 = d4^2 - d1^2 - a4^2 + a1^2 - b4^2 + b1^2 - c4^2 + c1^2
2ya1 - 2ya4 = d4^2 - d1^2 - a4^2 + a1^2 + b1^2 - b4^2 - c4^2 + c1^2
2za1 - 2za4 = d4^2 - d1^2 - a4^2 + a1^2 - b4^2 + b1^2 + c1^2 - c4^2

将上述方程表示为矩阵形式Ax = b，其中：

A = [2(a1-a2) 2(b1-b2) 2(c1-c2);
2(a1-a3) 2(b1-b3) 2(c1-c3);
2(a1-a4) 2(b1-b4) 2(c1-c4)];

x = [x; y; z];

b = [d2^2-d1^2-a2^2+a1^2-b2^2+b1^2-c2^2+c1^2;
d3^2-d1^2-a3^2+a1^2-b3^2+b1^2-c3^2+c1^2;
d4^2-d1^2-a4^2+a1^2-b4^2+b1^2-c4^2+c1^2];

解出x即可得到卡车的位置。

以下是MATLAB程序：

% 卫星位置（单位：千米）
a1 = 15.6; b1 = 7.54; c1 = 20.14;
a2 = 918.76; b2 = 2.75; c2 = 18.61;
a3 = 17.61; b3 = 14.63; c3 = 13.48;
a4 = 19.17; b4 = 0.61; c4 = 18.39;

% 卡车到卫星的距离（单位：千米）
d1 = 22.167124;
d2 = 22.604821;
d3 = 24.013845;
d4 = 22.670775;

% 构造矩阵A和向量b
A = [2*(a1-a2) 2*(b1-b2) 2*(c1-c2);
     2*(a1-a3) 2*(b1-b3) 2*(c1-c3);
     2*(a1-a4) 2*(b1-b4) 2*(c1-c4)];
b = [d2^2-d1^2-a2^2+a1^2-b2^2+b1^2-c2^2+c1^2;
     d3^2-d1^2-a3^2+a1^2-b3^2+b1^2-c3^2+c1^2;
     d4^2-d1^2-a4^2+a1^2-b4^2+b1^2-c4^2+c1^2];

% 解出x
x = A \ b;
x = x';

% 显示卡车的位置
fprintf('卡车的位置为：(%.2f, %.2f, %.2f)\n', x(1), x(2), x(3));

运行上述程序即可得到卡车的位置为：(1.88, 2.64, 0.46) 千米。