idnlhw/findop

Compute operating point for Hammerstein-Wiener model

Syntax

[X,U] = findop(sys,'steady',InputLevel,OutputLevel)

[X,U] = findop(sys,spec)

[X,U] = findop(___,Options)

[X,U,Report] = findop(___)

[X,U] = findop(sys,'snapshot',T,Uin)

[X,U] = findop(sys,'snapshot',T,Uin,X0)

Description

[X,U] = findop(sys,'steady',InputLevel,OutputLevel) returns the operating-point state values, X, and input values, U, for the idnlhw model, sys, using steady-state input and output specifications.

[X,U] = findop(sys,spec) returns the steady-state operating point for sys using the operating point specification in spec.

[X,U] = findop(___,Options) specifies optimization search options for all of the previous syntaxes.

[X,U,Report] = findop(___) returns a summary report on the optimization search results for all of the previous syntaxes.

[X,U] = findop(sys,'snapshot',T,Uin) returns the operating point for sys at a simulation snapshot at time, T, using the specified input, Uin. The initial states of sys are assumed to be zero.

[X,U] = findop(sys,'snapshot',T,Uin,X0) specifies the initial states of the simulation.

Examples

Find Steady-State Hammerstein-Wiener Operating Point Using Default Specifications

Load the estimation data and estimate a Hammerstein-Wiener model.

load twotankdata;
z = iddata(y,u,1);
M = nlhw(z,[5 1 3]);

Find the steady-state operating point where the input level is set to 1 and the output is unknown.

[X,U] = findop(M,'steady',1,NaN);

Find Hammerstein-Wiener Operating Point Using Additional Specifications

Estimate a Hammerstein-Wiener model.

load iddata7;
orders = [4*ones(1,2) 2*ones(1,2) 3*ones(1,2)];
M = nlhw(z7,orders,'unitgain','pwlinear');

Warning: "unitgain" is deprecated and may be removed in a future release. Use "idUnitGain" instead.

Warning: "pwlinear" is deprecated and may be removed in a future release. Use "idPiecewiseLinear" instead.

Create a default operating point specification object.

spec = operspec(M);

Set the values for the input signals.

spec.Input.Value(1) = -1;
spec.Input.Value(2) = 1;

Set the maximum and minimum values for the output signal.

spec.Output.Max = 10;
spec.Output.Min = -10;

Find the steady-state operating point using the given specifications.

[X,U] = findop(M,spec);

Find Hammerstein-Wiener Operating Point Using Custom Options

Load the estimation data and estimate a Hammerstein-Wiener model.

load twotankdata;
z = iddata(y,u,1);
M = nlhw(z,[5 1 3]);

Create a default findopOptions option set.

opt = findopOptions(M);

Modify the option set to specify a steepest descent gradient search method with a maximum of 50 iterations.

opt.SearchMethod = 'grad';
opt.SearchOptions.MaxIterations = 50;

Find the steady-state operating point using the specified options.

[X,U] = findop(M,'steady',1,NaN,opt);

Retrieve Hammerstein-Wiener Operating Point Search Report

Load the estimation data and estimate a Hammerstein-Wiener model.

load iddata7;
orders = [4*ones(1,2) 2*ones(1,2) 3*ones(1,2)];
M = nlhw(z7,orders,[],idPiecewiseLinear);

Find the steady-state operating point where input 1 is set to 1 and input 2 is unrestricted. The initial guess for the output value is 2.

[X,U,R] = findop(M,'steady',[1 NaN],2);

Display the summary report.

disp(R);

            SearchMethod: 'auto'
                 WhyStop: 'Near (local) minimum, (norm(g) < tol).'
              Iterations: 20
               FinalCost: 4.4373e-31
    FirstOrderOptimality: 2.4722e-16
            SignalLevels: [1x1 struct]

Find Hammerstein-Wiener Simulation Snapshot Using Default Initial States

Load the estimation data estimate a Hammerstein-Wiener model.

load twotankdata;
z = iddata(y,u,1);
M = nlhw(z,[5 1 3]);

Find the simulation snapshot after 10 seconds, assuming initial states of zero.

[X,U] = findop(M,'snapshot',10,z);

Find Hammerstein-Wiener Simulation Snapshot Using Initial State Specifications

Load the estimation data and estimate a Hammerstein-Wiener model.

load twotankdata
z = iddata(y,u,1);
M = nlhw(z,[5 1 3]);

Create an initial state vector.

X0 = [10;10;5;5;1;1;0];

Find the simulation snapshot after 10 seconds using the specified initial states.

[X,U] = findop(M,'snapshot',10,z,X0);

Input Arguments

sys — Hammerstein-Wiener model

idnlhw object

Hammerstein-Wiener model, specified as an idnlhw object.

InputLevel — Steady-state input level

vector

Steady-state input level for computing the operating point, specified as a vector. The length of InputLevel must equal the number of inputs specified in sys.

The optimization algorithm assumes that finite values in InputLevel are fixed input values. Use NaN to specify unknown input signals with initial guesses of 0. The minimum and maximum bounds for all inputs have default values of -Inf and +Inf respectively.

OutputLevel — Steady-state output level

vector

Steady-state output level for computing the operating point, specified as a vector. The length of OutputLevel must equal the number of outputs specified in sys.

The values in OutputLevel indicate initial guesses for the optimization algorithm. Use NaN to specify unknown output signals with initial guesses of 0. The minimum and maximum bounds for all outputs have default values of -Inf and +Inf respectively.

spec — Operating-point specifications

operspec object

Operating-point specifications, such as minimum and maximum input/output constraints and known inputs, specified as an operspec object.

T — Operating point snapshot time

positive scalar

Operating point snapshot time, specified as a positive scalar. The value of T must be in the range [T₀, N*T_s], where N is the number of input samples, T_s is the sample time and T₀ is the input start time (Uin.Tstart).

Uin — Snapshot simulation input

iddata object | matrix

Snapshot simulation input, specified as one of the following:

• Time-domain iddata object with a sample time and input size that matches sys.

• Matrix with as many columns as there are input channels. If the matrix has N rows, the input data is assumed to correspond to the time vector (1:N)*sys.Ts.

X0 — Initial states

column vector

Initial states of the simulation, specified as a column vector with length equal to the number of states in sys. X0 provides the initial conditions at the time corresponding to the first input sample (Uin.Start, if Uin is an iddata object, or sys.Ts if Uin is a double matrix).

For more information about the states of an idnlhw model, see Definition of idnlhw States.

Options — Operating point search options

findopOptions option set

Operating point search options, specified as a findopOptions option set.

Output Arguments

X — Operating point state values

column vector

Operating point state values, returned as a column vector of length equal to the number of model states.

U — Operating point input values

column vector

Operating point input values, returned as a column vector of length equal to the number of inputs.

Report — Search result summary

structure

Search result summary report, returned as a structure with the following fields:

Field	Description
SearchMethod	Search method used for iterative parameter estimation. See SearchMethod in findopOptions for more information.
WhyStop	Search algorithm termination condition.
Iterations	Number of estimation iterations performed.
FinalCost	Final value of the minimization objective function (sum of the squared errors).
FirstOrderOptimality	$$\infty$$-norm of the search gradient vector when the search algorithm terminates.
SignalLevels	Structure containing the fields Input and Output, which are the operating point input and output signal levels respectively.

Algorithms

findop computes the operating point from steady-state operating point specifications or at a simulation snapshot.

Computing the Operating Point from Steady-State Specifications

To compute the steady-state operating point, call findop using either of the following syntaxes:

[X,U] = findop(sys,'steady',InputLevel,OutputLevel)
[X,U] = findop(sys,spec)

findop uses a different approach to compute the steady-state operating point depending on how much information you provide for this computation:

• When you specify values for all input levels (no NaN values). For a given input level, U, the equilibrium state values are X = inv(I-A)*B*f(U), where [A,B,C,D] = ssdata(model.LinearModel), and f() is the input nonlinearity.

• When you specify known and unknown input levels. findop uses numerical optimization to minimize the norm of the error and compute the operating point. The total error is the union of contributions from e₁ and e₂ , e(t) = (e₁(t)e₂(t)), such that:

  • e₁ applies for known outputs and the algorithm minimizes e₁ = y- g(L(x,f(u))), where f is the input nonlinearity, L(x,u) is the linear model with states x, and g is the output nonlinearity.

  • e₂ applies for unknown outputs and the error is a measure of whether these outputs are within the specified minimum and maximum bounds. If a variable is within its specified bounds, the corresponding error is zero. Otherwise, the error is equal to the distance from the nearest bound. For example, if a free output variable has a value z and its minimum and maximum bounds are L and U, respectively, then the error is e₂= max[z-U, L-z, 0].

  The independent variables for the minimization problem are the unknown inputs. In the error definition e, both the input u and the states x are free variables. To get an error expression that contains only unknown inputs as free variables, the algorithm findop specifies the states as a function of inputs by imposing steady-state conditions: x = inv(I-A)*B*f(U), where A and B are state-space parameters corresponding to the linear model L(x,u). Thus, substituting x = inv(I-A)*B*f(U) into the error function results in an error expression that contains only unknown inputs as free variables computed by the optimization algorithm.

Computing the Operating Point at a Simulation Snapshot

When you use the syntax [X,U] = findop(sys,'snapshot',T,UIN,X0), the algorithm simulates the model output until the snapshot time, T. At the snapshot time, the algorithm computes the inputs for the linear model block of the Hammerstein-Wiener model (LinearModel property of theidnlhw object) by transforming the given inputs using the input nonlinearity: w = f(u). findop uses the resulting w to compute x until the snapshot time using the following equation: x(t+1) = Ax(t) + Bw(t), where [A,B,C,D] = ssdata(model.LinearModel).

Note

For snapshot-based computations, findop does not perform numerical optimization.

Extended Capabilities

Automatic Parallel Support

Accelerate code by automatically running computation in parallel using Parallel Computing Toolbox™.

Parallel computing support is available for estimation using the lsqnonlin search method (requires Optimization Toolbox™). To enable parallel computing, use findopOptions, set SearchMethod to 'lsqnonlin', and set SearchOptions.Advanced.UseParallel to true.

For example:

opt = findopOptions(idnlhw);
opt.SearchMethod = 'lsqnonlin';
opt.SearchOptions.Advanced.UseParallel = true;