How to properly compute weights for Weighted Least Squares (WLS)?

I want to apply the weighted least squares method in order to identify parameters of a dynamic process. The process is described by a second order differential equation of the form:

$$\ddot{y}+a_1\dot{y}+a_2y=b_0u$$

where $$\ y,u$$ are the output and input of the process respectively and are measured as well. The unknown parameters $$\ a_1,a_2,b_0$$ are those that have to be identified by processing the input-output data of the system. The process can also be written in the form:

$$y = [a_1 \ a_2 \ b_0]Ψ \rightarrow y=\theta^{*T}Ψ$$

where $$\ Ψ$$ is computed through the input-output data of the system.

I have already implemented the Ordinary Least Squares (OLS) method given by the formula:

$$\hat{\theta}= (Ψ^ΤΨ)^{-1}Ψ^Τy$$

so the matrix $$\ Ψ$$ is computed and $$\ \hat{\theta}=[\hat{a_1} \ \hat{a_2} \ \hat{b_0}]$$ are the estimates of the parameters. In order to implemented the weighted method I need to apply the formula:

$$\hat{\theta}= (Ψ^ΤWΨ)^{-1}Ψ^ΤWy$$

where $$\ W$$ is a diagonal matrix. I have searched some literature and found out that the "best" choice of $$\ W$$ is $$\ W = (E[ee^T])^{-1}$$ where $$\ e = y - \hat{y}$$ ($$\ \hat{y}$$ is the estimated output of the process). However, the covariance matrix of the error equation is not known and this form of $$\ W$$ can't be obtained and used. I also haven't found any good guidance (or an algorithmic one) in order to properly compute the weights' matrix. Is there an algorithmic way by using the vectors $$\ y,u$$ and the matrix $$\ Ψ$$ in order to properly compute $$\ W$$ and apply the WLS method ?

In your case, you would use that if the parameters $$y$$, $$\dot y$$, and $$\ddot y$$ deviated by several orders of magnitude. However, you would do that then before the least squares (which the induces the diagonal matrix $$W$$).