Simple optimizer?

I have the following equation:

res = c1 * v_app * u_s + c2/v_app * sin(psi) * cos(beta) - psi_dot

and I want to determine the constants c1 and c2 based on measurements such that the the sum of res² over all measurements reaches a minimum.

I have vectors of v_app, u_s, psi, beta and psi_dot with about 1000 elements.

Which package would you suggest?

Re-express as a linear system Ac = b where c=[c_1, c_2] is a vector of the unknown coefficients, A is a 1000 \times 2 matrix whose columns are v_{app} u_s and \sin(\psi) \cos(\beta) /v_{app}, and b is the vector psi_dot.

Then solve with c = A\b. That’ll get you a least-squares solution to the problem, the one that minimizes sum res^2, using QR decomposition.

5 Likes

Very nice idea! No package needed at all!

Works nicely:

function calc_c1_c2(v_app, psi, beta, psi_dot, steering)
    col1 = v_app .* steering
    col2 = sin.(psi) .* cos.(beta) ./ v_app
    A = [col1 col2]
    c = A \ psi_dot
    return c[1], c[2]
end
1 Like

Extra question: How can I obtain the standard deviation of c1 and c2 ?

You can use the StatsBase package’s std(c1) function.

Not really. If I calculate:

c = A \ psi_dot

I get a two element vector of c1 and c2, both scalars. And I cannot calculate the standard deviation of a scalar.

The package LsqFit.jl has a standard_error() of the fit parameters.
See the tutorial section here.

But isn’t LsqFit intended for non-linear models, and here we have a linear model?

Absolutely not, performing linear regressions with LsqFit.jl is perfectly fine, and it is probably when the estimated standard deviations are the most accurate.

I think we are actually talking about a standard error of a parameter estimate here, not a standard deviation. You first should think about conceptually what you are actually trying to estimate here, as the docs of LsqFit say there are some assumptions being made

Assume the errors in each sample are independent, normal distributed with zero mean and same variance, i.e. ϵ∼N(0,σ2I)

which might or might not make sense (especially homoskedasticity is often a problematic assumption on real world data at least in the social sciences, I’m not a kite power plant modelling expert so YMMV).

But indeed your model is linear so you can just estimate it with GLM or FixedEffectModels and get heteroskedasiticty-robust standard errors as well if that matters to you.

1 Like

Some terminology reference.

1 Like

what if you replace the terms in your A and psi_dot arrays for Measurements.measurement types?