Special TF Optimization in CAPTAIN

Friday, 29 April 2011

For special TF models, such as those constrained to have real eigenvalues (poles) or Hammerstein models, it is possible to embed the rivbj or rivcbj routines into a Matlab lsqnonlin routine. For example, in the former case, a typical optimization routine would be:

function Q=CTconstrainedopt(X,y,u,m,disc)

A=conv([1 X(1)],[1 X(2)]);A=conv(A,[1 X(3)]);% constrain denominator A to have real poles

if disc

 [th,M,D,stats,e,eef]=...

 rivbj([y u],m,[0.00001 -1 -1 1],...

 [A(2:end) 0 0 0],diag([0 0 0 10^6 10^6 10^6]));% recursive estimation with poles constrained to A above

else

 [th,M,Mn,D,stats,e,eef]=...

 rivcbj([y u],m,[0.00001 2 1 0],...

 -0.1,[A(2:end) 0 0 0],diag([0 0 0 10^6 10^6 10^6]));% recursive estimation with poles constrained to A above

end

 xh=y-e;plot([y xh]);drawnow; % not essential just shows how model is explaining the data

 Q=e;% criterion (could be eef if noise model is being estimated concurrently)

 if ~isfinite(Q);Q=1e+23;end

Here, the recursive option of the rivbj or rivcbj routines is selected and the initial estimates of the TF denominator parameters are set to the latest optimized values of the A denominator coefficients, with the associated initial diagonal elements of the covariance matrix set to zero in order to constrain them to these latest optimized values. The initial estimates of the numerator B parameters are set to zero and the associated initial diagonal elements of the covariance matrix set to 10^6 in order to allow them to be estimated in the normal manner. The covariance matrix of the B parameter estimates is recovered from the rivbj or 

rivcbj routine; and the covariance matrix of the constrained poles (here X(1), X(2) and X(3)) is obtained by calling the CAPTAIN confint routine. 

In the case of the Hammerstein model, the optimized parameter vector X would contain the parameters of the input nonlinearity and the input u wouldbe defined by the nonlinear function using these latest optimized estimates as ue. The rivbj and rivcbj estimation would then be unconstrained with the input variable defined by the latest nonlinearly transformed input ue.