How do I minimize a particular error function with lsqnonlin?

3 ビュー (過去 30 日間)
Debdeeprc
Debdeeprc 2021 年 2 月 23 日
コメント済み: Debdeeprc 2021 年 2 月 24 日
Hello all,
I am trying to minimize the error function
E = abs(vecnorm(b_hat) - vecnorm(b_ref))
by solving for A and x in
b_hat = A * (b_raw - d)
where b_ref, b_raw and b_hat have dimensions of 3 x 100 (they are all 3 dimensional vectors, with 100 different values over time). A is a 3 x 3 matrix and d is a 3 x 1 vector. Only the values of b_raw and vecnorm(b_ref) are known to me. The correct values of A and d would help me to calculate b_hat, whose norm would be equal to that of b_ref. At the moment, I am solving them using lsqnonlin with the option 'trust-region-reflective', which helps me also to set up the upper and lower bounds, like this:
fun = @(x) abs(vecnorm([x(1) x(2) x(3); x(4) x(5) x(6); x(7) x(8) x(9)] * (b_raw - [x(10); x(11); x(12)])) - vecnorm(b_ref));
x0(1, 1:9) = [1 0 0 0 1 0 0 0 1];
x0(1, 10:12) = [0 0 0];
options = optimoptions('lsqnonlin');
options.Algorithm = 'trust-region-reflective';
lb = [0.8 -0.1 -0.1 -0.1 0.8 -0.1 -0.1 -0.1 0.8 -100000 -100000 -100000];
ub = [1.2 0.1 0.1 0.1 1.2 0.1 0.1 0.1 1.2 100000 100000 100000];
options.StepTolerance = 1e-10;
options.FunctionTolerance = 1e-10;
options.OptimalityTolerance = 1e-10;
options.MaxFunctionEvaluations = 1e10;
options.MaxIterations = 1e3;
[y, resnorm] = lsqnonlin(fun, x0, lb, ub, options);
A = [y(1) y(2) y(3);...
y(4) y(5) y(6);...
y(7) y(8) y(9)];
d = [y(10); y(11); y(12)];
Is there a better way to solve this? If I set the limits far-apart too much, the estimation of A and d becomes incorrect.
Thanks a lot!
  3 件のコメント
Debdeeprc
Debdeeprc 2021 年 2 月 23 日
Hi Matt,
Thanks for your comment. Please note, A * b_raw-d is not even close to b_ref. Only, their norms are "almost equal". I will give more context here to clarify. b_raw is actually the magnetic field output vector from an uncalibrated magnetic field sensor which is on a low earth orbiting satellite. The sensor output values (b_raw) have to be corrected using a 3 x 3 calibration matrix A and an offset vector d to get the calibrated output b_hat. The information about what the correct magnetic field norm should be is already available from the IGRF magnetic field model, which in my case depends on the satellite's position. The IGRF gives the output in an Earth Centred frame (b_ref). A calibrated magnetic field sensor output should have the same norm/magnitude as b_ref, however the individual vector elements would be different as they are in two different reference frames: b_hat is in satellite body fixed frame and b_ref is Earth centred.
Thanks again!
Debdeeprc
Debdeeprc 2021 年 2 月 23 日
Hi Matt,
A is NOT a rotation matrix. It is more commonly referred to as a calibration matrix, whose diagonal terms scale the raw sensor values and the off-diagonal terms correct the orthogonality errors in the magnetic field sensor (sometimes cheaper sensors don't have their axes mutually perpendicular to one another). Actually, you can read more about this specific part under the section Correction Technique of https://www.mathworks.com/help/nav/ug/magnetometer-calibration.html
Thanks!

サインインしてコメントする。

採用された回答

Matt J
Matt J 2021 年 2 月 23 日
編集済み: Matt J 2021 年 2 月 24 日
I don't see much that I think can be improved, except perhaps to formulate it so that your errors are differentiable - the lsqnonlin algorithms are generally driven by Jacobian calculations.
v_ref=vecnorm(b_ref).^2; %compute only once
B=[b_raw; -ones(1,100)]; %put into homogeneous form
fun = @(Ad) vecnorm(Ad*B).^2 - v_ref; %Use difference of vecnorm^2 for differentiability
options = optimoptions('lsqnonlin');
options.Algorithm = 'trust-region-reflective';
lb = [0.8 -0.1 -0.1 -0.1 0.8 -0.1 -0.1 -0.1 0.8 -100000 -100000 -100000];
ub = [1.2 0.1 0.1 0.1 1.2 0.1 0.1 0.1 1.2 100000 100000 100000];
options.StepTolerance = 1e-10;
options.FunctionTolerance = 1e-10;
options.OptimalityTolerance = 1e-10;
options.MaxFunctionEvaluations = 1e10;
options.MaxIterations = 1e3;
[Ad, resnorm] = lsqnonlin(fun, eye(3,4), lb, ub, options);
A=Ad(1:3,1:3);
d=Ad(:,4);
  1 件のコメント
Debdeeprc
Debdeeprc 2021 年 2 月 24 日
Hi Matt,
Thanks a lot for your suggestions. I already see improvements in the results after implementing your suggestions, especially in cases where sufficient variation in vecnorm(b_ref) is not available. Using the difference of vecnorm^2 instead of vecnorm works better.
Have a great day,
Debdeep

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeSatellite and Orbital Mechanics についてさらに検索

製品


リリース

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by