Advice on formulating a minimization function

4 次查看（过去 30 天）
Collin Pecora
Collin Pecora 2024-2-1
编辑： William Rose 2024-2-3
I am modeling ballistic data using the Low-Angle Trajectory approximation of quadratic drag from
R. D. H. Warburton and J. Wang, “Analysis of asymptotic projectile motion with air resistance using the Lambert W function,” American Journal of Physics, Vol. 72,pp. 1404–1407, 2004.
the revelant functional forms of (5) and (6) are
vxFcn = @(v0,th,b,t) v0*cos(th)./(1+v0*cos(th)*b.*t);
vyFcn = @(v0,th,b,t) ((v0*sin(th)-0.5*g.*t)./(1+v0*cos(th)*b.*t)) - 0.5*g.*t;
xFcn = @(v0,th,b,t) (1/b) * log(1 + v0*b*cos(th)*t);
yFcn = @(v0,th,b,t) (v0*sin(th) + (g/(2*b*v0*cos(th)))).*(1/(b*v0*cos(th))).*log(1+b*v0*cos(th).*t)-(0.25*g.*t.^2);
The measurements of vx and vy are noisy, while the measured impact locations, ximpact and yimpact are precise.
I use fmincon twice, the first pass uses velocityMinFcn below, with
k = [v0,th,b]. Here t is the vector of measurement times
function SSE = velocityMinFcn(k)
vxm = vxFcn(k(1),k(2),k(3),t);
vym = vyFcn(k(1),k(2),k(3),t);
SSE = sum((vxm-vx).^2) + sum((vym-vy).^2);
end
I then do a refinment pass, again using fmincon, with k = [v0,th,b,impactTime] using v0, th and b from the first pass
function SSE = impactLocationFcn(k)
SSE = abs(xFcn(k(1),k(2),k(3),k(4)) - ximpact);
SSE = SSE + abs(yFcn(k(1),k(2),k(3),k(4)) - yimpact);
end
My question is I would like to combine the two minimization function into one and am looking for advice on how best to do it
  2 个评论
William Rose
William Rose 2024-2-1
编辑：William Rose 2024-2-1
@Collin Pecora,
Please provide some example data, and the script you are using thus far, if you have one.
Am I right to understand that you are trying to esitmate v0, theta, and b (drag coiefficient) that is consistent with the measured data?
Am I right to understand that you have measured data:
  • vectors x, y, vx, vy, t
  • scalars ximpact, yimpact, timpact
Are the functions vxFcn, vyFcn, xFcn, yFcn supposed to return vectors when called with vector time, and scalar when called with scalar time?
When you call impactLocationFcn, do you pass a scalar time for k(4)?
Your function SSE=impactLocationFcn refers to ximpact and to yimpacty. yimpacty looks like a spelling error which could cause an error.
Collin Pecora
Collin Pecora 2024-2-1
Thanks for the reply
I have the measured, noisy, data , vx, yy and t which are vectors
I have the measured scalars, ximpact, yimpact
k(4), impact time, is a parameter to be solved by fmincon in the refinement pass.
In the end, I am after v0, th, b and impact time.
Typo has been fixed.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

William Rose
William Rose 2024-2-1
@Collin Pecora,
If you want to use multiple sources of data, with different levels of measurement error, to estimate parameters, then weight the errors by a factor proportional to , where σ is the standard deviation of the measurement for that quantity. In other words, suppose you think range (R=final distance) is 5 times more accurate than velocity (vx, vy, in the units of measure: since position and velocity have different units, their respective sigmas shoudl be measured in the respecitve units). Then the weighing factors for range and velocity are .
The weighted sum squared error, which you want to minimize, is
where subscripts e and m are for estimated and measured, and there are N measurements of velocity.
Call fmincon() once, to estimate the value of [v0,theta,b] that minimizes wSSE.
If you are not sure about the sigmas, then make your best guess.
I would not use fmincon to solve for the impact time. Based on your comment, you do not have a measured impact time to compare to a predicted impact time. Compute the impact time using the estimates of v0, theta,b which you will get from fmincon(). I assume the elevation (y-value) of the launch and the landing spot are the same, or if not the same, the difference is known. If the landing zone is a hill, the landing elevation will be range-dependent, which complicates things.
  3 个评论
显示 1更早的评论
Collin Pecora
Collin Pecora 2024-2-2
移动：Torsten 2024-2-2
@William Rose
I was lazy and didn't update the reference, they are (5) and (6) from
R. Warburton, J. Wang and J. Burgdörfer, "Analytic Approximations of Projectile Motion with Quadratic Air Resistance," Journal of Service Science and Management, Vol. 3 No. 1, 2010, pp. 98-105. doi: 10.4236/jssm.2010.31012.
Collin
William Rose
William Rose 2024-2-2
编辑：William Rose 2024-2-3
@Collin Pecora,
[Edit: Upload revised script projectileMotionFit. New version is shorter, clearer, gves same results.]
I used the W&W 2004 equations, before I noticed your post about W&W 2010. I simulated the system and added Gaussian random noise to the positions (which includes range) and to the velocities, . I used the same initial conditions as W&W 2004, Figure 2, and my results look like their Fig 2, but my results are noisy, unlike theirs. See below.
I saved the results in file data1.mat, attached. The file has 15 variables: t1,x1,y1,vx1,vy1 (red data above, ); t2,...,vy2 (green data, ); and t3,...,vy3 (blue data,).
Script projectileMotionFit.m reads the noisy data from the data file. It determines the range from the position data. It estimates [v0,theta,b] by minimizing the squared error. The squared error function uses sigmas for inverse weighting: sigma=1 for range and sigma=10 for velocities. The true [v0,theta,b] are [300,30,1.00]. Results below:
>> projectileMotionFit
Range 1=260.1, R2=281.5, R3=294.9 m.
Fit 1: v0=298.5 m/s, theta=30.4 deg, b=0.991 /s.
Fit 1: Vel.rmse=10.0 m/s, rangeErr=-0.3 m.
>>
W&W 2004 say "This paper was inspired by the results of an honors introductory course", referencing Mr. David Pierce. That must have been a later-year-version of the high school AP physics course I took with him. Thank you, Mr. Pierce!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Curve Fitting Toolbox 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by