利用c#实现磁偶极子定位时，采用梯度下降法，优化后的磁偶极子的位置和磁距竟然与设置的初始值是一样的，但是不对啊，优化的结果应该跟初值没关系才对，附完整代码，求助问题出在哪

我观此物与我有缘 2024-11-25 16:19:38

using System;

//基于梯度下降法的磁偶极子定位算法——完整

class Complete_Magnetic_Dipole_Localization_Algorithm
{
// 真空磁导率
const double Mu0 = 4 * Math.PI * 1e-7;
// 假设观测到的磁场——已知

// 计算磁场

static double[][] ComputeMagneticField(double[][] sensorPositions, double[] dipolePosition, double[] dipoleMoment)
{
int numSensors = sensorPositions.Length;
double[][] magneticFields = new double[numSensors][];
for (int i = 0; i < numSensors; i++)
{
// 获取传感器位置
double[] sensorPosition = sensorPositions[i];
double[] r = new double[3];
for (int j = 0; j < 3; j++)
{
r[j] = sensorPosition[j] - dipolePosition[j];
}
//计算r的平方和
double rSquared = 0;
for (int j = 0; j < 3; j++)
{
rSquared += r[j] * r[j];
}
//计算r的三次方
double rCubed = Math.Pow(rSquared, 1.5);
//计算磁距与向量的点积
double dotProduct = 0;

for (int j = 0; j < 3; j++)
{
dotProduct += dipoleMoment[j] * r[j]; // 计算磁矩与位置向量的点积
}
//计算磁场
double[] magneticField = new double[3];
for (int j = 0; j < 3; j++)
{
magneticField[j] = (Mu0 / (4 * Math.PI)) * ((3 * dotProduct * r[j]) / rCubed - dipoleMoment[j] / (rSquared * Math.Sqrt(rSquared)));
}
// 存储每个传感器位置的磁场
magneticFields[i] = magneticField;
}

return magneticFields; //计算B
}

// 计算磁场与观测磁场之间的误差（损失函数）
static double ComputeLoss(double[][] sensorPositions, double[][] observedFields, double[] initialdipolePosition, double[] initialdipoleMoment)
{
double loss = 0;
for (int i = 0; i < sensorPositions.Length; i++)//for (int i = 0; i < sensorPositions.Length; i++)
{
double[] predictedField = ComputeMagneticField(new double[][] { sensorPositions[i] }, initialdipolePosition, initialdipoleMoment)[0];
for (int j = 0; j < 3; j++)
{

loss += Math.Pow(predictedField[j] - observedFields[i][j], 2);
}

}
return loss;
}

//坐标变换
public static double[][] TransformCoordinates(double[][] positions, int originIndex)
{
double[][] transformedPositions = new double[positions.Length][];
double[] origin = positions[originIndex]; // 选择某个传感器作为原点

for (int i = 0; i < positions.Length; i++)
{
transformedPositions[i] = new double[3];
for (int j = 0; j < 3; j++)
{
transformedPositions[i][j] = positions[i][j] - origin[j]; // 相对原点进行坐标变换
}
}

return transformedPositions;
}

// 计算磁偶极子位置的梯度
static double[] ComputeGradientWithRespectToPosition(double[][] sensorPositions, double[][] observedFields, double[] initialdipolePosition, double[] initialdipoleMoment)
{
double[] gradient = new double[3];
double epsilon = 1e-6;
for (int i = 0; i < 3; i++)
{
double[] perturbedPosition = (double[])initialdipolePosition.Clone(); // 复制当前的位置
perturbedPosition[i] += epsilon; // 在第i个方向上加上一个小扰动（epsilon）

double lossPlus = ComputeLoss(sensorPositions, observedFields, perturbedPosition, initialdipoleMoment);

double lossMinus = ComputeLoss(sensorPositions, observedFields, initialdipolePosition, initialdipoleMoment);

gradient[i] = (lossPlus - lossMinus) / epsilon;
}
return gradient;
}

// 计算磁偶极矩的梯度
static double[] ComputeGradientWithRespectToMoment(double[][] sensorPositions, double[][] observedFields, double[] initialdipolePosition, double[] initialdipoleMoment)
{
double[] gradient = new double[3];
double epsilon = 1e-6;

for (int i = 0; i < 3; i++)
{
double[] perturbedMoment = (double[])initialdipoleMoment.Clone();
perturbedMoment[i] += epsilon;

double lossPlus = ComputeLoss(sensorPositions, observedFields, initialdipolePosition, perturbedMoment);
double lossMinus = ComputeLoss(sensorPositions, observedFields, initialdipolePosition, initialdipoleMoment);

gradient[i] = (lossPlus - lossMinus) / epsilon;
}
return gradient;
}

// 梯度下降算法优化磁偶极子的位置和磁矩
static (double[], double[]) GradientDescentOptimization(double[][] sensorPositions, double[] initialdipolePosition, double[] initialdipoleMoment, double[][] observedFields, double learningRate, int iterations)
{
//double tolerance = 1e-6;
for (int iter = 0; iter < iterations; iter++)
{
// 计算损失
double loss = ComputeLoss(sensorPositions, observedFields, initialdipolePosition, initialdipoleMoment);

double[] positionGradient = ComputeGradientWithRespectToPosition(sensorPositions, observedFields, initialdipolePosition, initialdipoleMoment);
double[] momentGradient = ComputeGradientWithRespectToMoment(sensorPositions, observedFields, initialdipolePosition, initialdipoleMoment);

double decay = Math.Pow(0.9, iter / 100); // 每100次迭代减小10%
double initialLearningRate = learningRate;
double adjustedLearningRate = initialLearningRate * decay;
for (int i = 0; i < 3; i++)
{
initialdipolePosition[i] -= adjustedLearningRate * positionGradient[i];
initialdipoleMoment[i] -= adjustedLearningRate * momentGradient[i];
}
// 每100次输出一次损失值
if (iter % 100 == 0)
{
Console.WriteLine($"Iteration {iter}, Loss: {loss}");
}
}
return (initialdipolePosition, initialdipoleMoment);
}

public static void Main(string[] args)
{

double[][] sensorPositions =
{

new double[] { 1.0, 2.0, 3.0 },
new double[] { 2.0, 3.0, 4.0 },
new double[] { 3.0, 4.0, 5.0 },
};

double[] dipolePosition = new double[] { 2.0, 2.0, 2.0 };
double[] dipoleMoment = new double[] { 2.0, 2.0, 3.0 };

double[][] observedFields = ComputeMagneticField(sensorPositions, dipolePosition, dipoleMoment);

double learningRate = 1e-4;
int iterations = 1000;//迭代次数超过40会导致损失值变大
double[] initialdipolePosition = new double[] { 1.0, 1.0, 1.0 };
double[] initialdipoleMoment = new double[] { 1.0, 1.0, 1.0 };

var (optimizedPosition, optimizedMoment) = GradientDescentOptimization(
sensorPositions, initialdipolePosition, initialdipoleMoment, observedFields, learningRate, iterations);

Console.WriteLine("Optimized Dipole Position: ");
Console.WriteLine($"X: {optimizedPosition[0]}, Y: {optimizedPosition[1]}, Z: {optimizedPosition[2]}");

Console.WriteLine("Optimized Dipole Moment: ");
Console.WriteLine($"X: {optimizedMoment[0]}, Y: {optimizedMoment[1]}, Z: {optimizedMoment[2]}");

}

...全文