I managed to get the method to work in JuMP.

Here is my code for the first example of Lars Nielsen’s paper:

```
using NamedArrays
using JuMP
import Ipopt
function lsq()
model = Model(Ipopt.Optimizer)
# measured input quantities
scale_indications = [ # [div]
199.988617
199.988608
174.992133
175.009992
150.013558
149.980675
125.002083
124.984217
100.005650
99.982925
74.986433
75.004325
50.007892
49.974992
24.978533
24.996417
]
scale_indications_reference = [199.998867, 199.998875] # [div]
weights_sum = 199.988816 # [g]
weight_reference = 199.999043 # [g]
density = 7965.76 # [kg/m³]
density_reference = 7833.01 # [kg/m³]
density_air = 1.1950 # [kg/m³]
# model input variables
@variables(model, begin
Iₒ[i=1:16] == scale_indications[i]
Irₒ[i=1:2] == scale_indications_reference[i]
msₒ == weights_sum
mrₒ == weight_reference
ρₒ == density
ρrₒ == density_reference
aₒ == density_air
# ζ parameters (refined from z)
I[i=1:16], (start = scale_indications[i])
Ir[i=1:2], (start = scale_indications_reference[i])
ms, (start = weights_sum)
mr, (start = weight_reference)
ρ, (start = density)
ρr, (start = density_reference)
a, (start = density_air)
# β parameters (unknown)
f, (start = 1.0) # scale correction factor
A, (start = 1.0) # non-linearity
m1, (start = 1.0) # mass 1 to 4
m2, (start = 1.0)
m3, (start = 1.0)
m4, (start = 1.0)
end
)
# collect input variables in array
zs = [[Iₒ[i] for i in 1:16]..., Irₒ[1], Irₒ[2], msₒ, mrₒ, ρₒ, ρrₒ, aₒ]
ζs = [ [I[i] for i in 1:16]..., Ir[1], Ir[2], ms, mr, ρ, ρr, a ]
# variance of input parameters
var_z = [5.29e-10*ones(18)..., 1e-10, 6.4e-11, 0.0841, 0.5041, 1.225e-5]
nz = length(zs)
# covariance matrix
# (named array allows expressive indexing)
Σ = NamedArray(zeros((nz, nz)), (zs, zs))
for (zi, var) in zip(zs, var_z)
Σ[zi, zi] = var
end
invΣ = inv(Σ) # fortunately, inversion leaves row and column names alone
# objective
# minimize difference of measured and estimated parameters
# weighted according to (co)variance
@NLobjective(model, Min, sum(
(zₗ - ζₗ)*invΣ[zₗ, zₖ]*(zₖ - ζₖ)
for (zₗ, ζₗ) in zip(zs, ζs)
for (zₖ, ζₖ) in zip(zs, ζs)
if invΣ[zₗ, zₖ] != 0
)
)
# constraints
mass_combinations = [
m1 + m2 + m3 + m4 # I1
m1 + m2 + m3 + m4 # I2
m1 + m2 + m3 # I3
m1 + m2 + m4 # I4
m1 + m2 # I5
m1 + m3 + m4 # I6
m1 + m4 # I7
m1 + m3 # I8
m1 # I9
m2 + m3 + m4 # I10
m2 + m3 # I11
m2 + m4 # I12
m2 # I13
m3 + m4 # I14
m3 # I15
m4 # I16
]
for i in 1:16
@NLconstraint(model, mass_combinations[i]*(1 - (a - 1.2)*(1/ρ - 1/8000)) - f*(I[i] + A*I[i]^2) == 0)
end
@NLconstraint(model, mr*(1 - (a - 1.2)*(1/ρ - 1/8000)) - f*(Ir[1] + A*Ir[1]^2) == 0)
@NLconstraint(model, mr*(1 - (a - 1.2)*(1/ρ - 1/8000)) - f*(Ir[2] + A*Ir[2]^2) == 0)
@NLconstraint(model, ms - (m1 + m2 + m3 + m4) == 0)
# abracadabra!
optimize!(model)
return(model)
end
mod = lsq()
solution_summary(mod; verbose=true)
```

The summary of the solution is:

```
* Solver : Ipopt
* Status
Termination status : LOCALLY_SOLVED
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Result count : 1
Has duals : true
Message from the solver:
"Solve_Succeeded"
* Candidate solution
Objective value : 8.336278266129565
Primal solution :
A : -4.29028123993008e-9
I[10] : 99.98289827591923
I[11] : 74.9864495677386
I[12] : 75.00432481389886
I[13] : 50.00788146321463
I[14] : 49.974995368628
I[15] : 24.978557386324223
I[16] : 24.996432624814283
I[1] : 199.98861731519125
I[2] : 199.98861731519125
I[3] : 174.9921471573928
I[4] : 175.0100224188919
I[5] : 150.0135576186036
I[6] : 149.98067149579722
I[7] : 125.0020873177321
I[8] : 124.98421206390323
I[9] : 100.00563324334139
Ir[1] : 199.99885452781737
Ir[2] : 199.99885452781737
Irₒ[1] : 199.998867
Irₒ[2] : 199.998875
Iₒ[10] : 99.982925
Iₒ[11] : 74.986433
Iₒ[12] : 75.004325
Iₒ[13] : 50.007892
Iₒ[14] : 49.974992
Iₒ[15] : 24.978533
Iₒ[16] : 24.996417
Iₒ[1] : 199.988617
Iₒ[2] : 199.988608
Iₒ[3] : 174.992133
Iₒ[4] : 175.009992
Iₒ[5] : 150.013558
Iₒ[6] : 149.980675
Iₒ[7] : 125.002083
Iₒ[8] : 124.984217
Iₒ[9] : 100.00565
a : 1.195
aₒ : 1.195
f : 1.0000018230343597
m1 : 100.00577238065725
m2 : 50.00796176585021
m3 : 24.978600179154746
m4 : 24.996475446351376
mr : 199.99904698570728
mrₒ : 199.999043
ms : 199.98880977201358
msₒ : 199.988816
ρ : 7965.76
ρr : 7833.01
ρrₒ : 7833.01
ρₒ : 7965.76
Dual solution :
* Work counters
Solve time (sec) : 0.00295
```

The results agree with the paper within the stated uncertainty.

My problem now is that I cannot figure out how to get the hessian at the solution, so I cannot calculate uncertainties and compare those.

Could anyone who is familiar with JuMP help me out here?