- vectors x, y, vx, vy, t
- scalars ximpact, yimpact, timpact
Advice on formulating a minimization function
4 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Collin Pecora
il 1 Feb 2024
Modificato: William Rose
il 3 Feb 2024
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 Commenti
William Rose
il 1 Feb 2024
Modificato: William Rose
il 1 Feb 2024
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:
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.
Risposta accettata
William Rose
il 1 Feb 2024
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 Commenti
William Rose
il 2 Feb 2024
Modificato: William Rose
il 3 Feb 2024
[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!
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Curve Fitting Toolbox in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!