> 
restart:
Digits:= 15:
St:= Statistics
:

> 
Zscore:= proc(V::Vector)
uses St=Statistics;
local mu:= St:Mean(V), sigma:= St:StandardDeviation(V);
mu, sigma, (V ~ mu)/sigma
end proc
:

> 
#Dependent variable:
Hvals:= <210.84, 218.19, 219.17, 222.27, 226.35, 228.45, 230.44, 234.65, 239.75>:
H__avg, H__sd, H__z:= Zscore(Hvals)
:

> 
#Independent variable:
Qvals:= <1800, 1750, 1730, 1700, 1650, 1635, 1600, 1570, 1500>:
Q__avg, Q__sd, Q__z:= Zscore(Qvals)
:

> 
#Proposed model, with (apparently) A known beforehand and Q__0, k__f, k__1, and Q__s to be
#determined.
model:=
subs(
Q= Q__N*Q__sigma + Q__mu,
H__N = (A*(1  Q/Q__0)  k__f*Q^2  k__1*(Q  Q__s)^2  H__mu)/H__sigma
)
:

> 
pv:= St:LinearFit(
[1,Q,Q^2], Q__z, H__z, Q, summarize,
output= [parametervector, residualsumofsquares]
);

Summary:

Model: .52025574e11.0044218*Q.58528770e1*Q^2

Coefficients:
Estimate Std. Error tvalue P(>t)
Parameter 1 0.0520 0.0414 1.2565 0.2556
Parameter 2 1.0044 0.0319 31.5038 0.0000
Parameter 3 0.0585 0.0325 1.8025 0.1215

Rsquared: 0.9941, Adjusted Rsquared: 0.9921
Solve for the symbolic paarameters in the model by equating coefficients of the powers of Q. We have 3 coefficients and 4 parameters, so I'll just set the value of k__1 (arbitrary choice) to the MathCAD value, 1e5.
> 
Known_params:=
[A, k__1, Q__mu, Q__sigma, H__mu, H__sigma]
=~ [350.156, 1e5, Q__avg, Q__sd, H__avg, H__sd]
;
PV:= solve(
eval(
{seq(pv[1]=~ coeff~(expand(rhs(model)), Q__N, [0,1,2]))},
Known_params
)
);

Compute the sum of the squared residuals. It should be the same as returned by the LinearFit command, allowing for some round off error.
> 
add(
eval(
eval(H__N  rhs(model), [PV[2][], Known_params[]])^2,
[H__N= H__z[k], Q__N= Q__z[k]]
),
k= 1..n
);

Compute the sum of the squared residuals using MathCAD's parameter values.
> 
add(
eval(
eval(
H__N  rhs(model),
[
([Q__0, k__f, Q__s]=~ [6962.32773, 2e5, 2198.86103])[],
Known_params[]
]
)^2,
[H__N= H__z[k], Q__N= Q__z[k]]
),
k= 1..n
);

