- vectors x, y, vx, vy, t
- scalars ximpact, yimpact, timpact

# Advice on formulating a minimization function

2 views (last 30 days)

Show older comments

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 Comments

William Rose
on 1 Feb 2024

Edited: William Rose
on 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.

### Accepted Answer

William Rose
on 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 Comments

William Rose
on 2 Feb 2024

Edited: William Rose
on 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!

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!