Advice on formulating a minimization function

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

[William Rose]

@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]

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.

Answer 1

[William Rose]

@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.

[William Rose]

@Collin Pecora,

I looked up the paper you cited, Warburton & Wang 2004. Warburton & Wang thank David Pierce for bringing the problem to their attention. He was my high school math and physics teacher, and an excellent teacher.

Your equations vxFcn and vyFcn do not match equations 2 and 4 of W&W 2004. Why not? What am I missing?

[Collin Pecora]

@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]

@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!