Getting Started

This section is created by exporting from the Matlab Live Script document docs/mat_source/gettingStarted.mlx

Feel free to download and work through that file in MATLAB or just read it here.

Getting started with iqcToolbox

Welcome! After this script, you will know how to install iqcToolbox, understand some of its basic commands, and hopefully have gotten a glimpse of how powerful these analysis tools are.

Table of Contents

Download Web browser MATLAB websave MATLAB Add-Ons Git Installation Building and analyzing uncertain systems Analyzing a system with no uncertainties Analyzing an uncertain system: Part I Analyzing an uncertain system: Part II Analyzing an uncertain system: Part III Reachability Analysis Incorporating Disturbances Incorporating uncertain initial conditions Conclusion References

Download

This toolbox can be downloaded a variety of ways:

Web browser

Enter https://github.com/iqcToolbox/iqcToolbox/archive/refs/heads/master.zip in your browser and unzip the downloaded package.

MATLAB websave

mkdir('path_to_iqcToolbox')

websave('iqcToolbox.zip', 'https://github.com/iqcToolbox/iqcToolbox/archive/refs/heads/master.zip')

unzip('iqcToolbox.zip', 'path_to_iqcToolbox')

MATLAB Add-Ons

In the main MATLAB window, go to Add-Ons -> Get Add-Ons and search for "iqcToolbox" and add it to MATLAB

Git

>> git clone https://github.com/iqcToolbox/iqcToolbox.git

NOTE:

Since this package is intended to continue developing open-source, users may prefer to git clone, or rely on the MATLAB add-on rather than download static archive files. Unfortunately, the MATLAB add-on portal does not update immediately after the GitHub repository is updated; there is typically a day lag between the MATLAB add-on version and the latest commit on GitHub. If you desire rapid updates, it is advisable to git clone the repository. Doing so will make updating this package much easier, via typical git pull/push commands.

Installation

This toolbox is installed with two easy commands:

addpath(fullfile('path_to_iqcToolbox', 'scripts'))

installIqcToolbox

installIqcToolbox will both download the necessary dependencies, and call initializeIqcToolbox, which will add the paths to those dependencies and the directories in iqcToolbox.

NOTE:

These two steps are necessary even if you obtained this package as a MATLAB add-on. Without ensuring certain dependencies are installed and added to the path, this toolbox will not correctly function.

NOTE:

Unless otherwise specified, installIqcToolbox will assume that:

• you want to interactively approve of which dependencies are downloaded
• you do NOT want to your added paths automatically saved for future MATLAB sessions; continuing with this default choice means you will need to call initializeIqcToolbox at the start of every MATLAB session

If you want to either of these options changed, call installIqcToolbox with both input arguments:

save_added_paths = true;

interactive = false;

installIqcToolbox(save_added_paths, interactive)

NOTE:

If you've already installed iqcToolbox, but now want to have all the correct paths permanently saved:

clear all

cfg_name = which('install_configuration.mat')

load(cfg_name)

save_added_paths = true;

save(cfg_name)

initializeIqcToolbox

Depending on your system, you may need admin rights for the added paths to be saved. For example, in a Linux system, start a matlab session with

>> sudo matlab

and repeat the previous MATLAB commands.

NOTE:

The only necessary dependencies for this toolbox are yalmip, SDPT3, and LPSOLVE. You may benefit from having other Semidefinite Program (SDP) solvers besides SDPT3. This might include SeDuMi, mosek, or csdp. The latter two solvers are capable of leveraging parallel computing resources; mosek is free for academic use, csdp is free for general use (EPL 2.0 license).

Building and analyzing uncertain systems

Let's start using iqcToolbox's most important functions.

Analyzing a system with no uncertainties

Consider the system :

Since the A-matrix is Schur (all eigenvalues are within the unit circle), we know the system is stable. Let's check this:

a = [0, 3/4; 3/4, 0];

b = [1; 0];

c = [1, 0];

d = 0;

timestep = -1; % An unspecified timestep

g = ss(a, b, c, d, timestep);

stable = isstable(g)

performance = norm(g, 'inf') % Measures H-infinity norm, or l2-induced performance, of stable sys

We can get this same information by converting the ss object g into a Ulft object and conducting IQC analysis on it:

g_lft = toLft(g)

% or

g_lft = toLft(a, b, c, d, timestep)

result = iqcAnalysis(g_lft)

result.valid % This boolean expresses if the LFT is robustly stable

result.performance % Returns upper-bound on a LFT's specified performance (default: l2-induced)

Analyzing an uncertain system: Part I

But what if there was uncertainty in the system model? Consider the case where our A-matrix could have some error in it's values:

We could prove that this is stable for all possible values of delta (robustly stable) by solving for the eigenvalues of the A-matrix and showing they are within the unit-circle. IQC analysis can also prove robust stability with the following:

dims = 1;

lower_bound = -1/2;

upper_bound = 1/2;

% a 1x1 Static, Linear, Time-Invariant (SLTI) uncertainty with the given bounds

del = DeltaSlti('del', dims, lower_bound, upper_bound)

a_del = [0, 3/4 + del; 3/4 - del, 0]; % the uncertain A-matrix

g_del = toLft(a_del, b, c, d, timestep) % the uncertain system defined with the uncertain A-matrix

result = iqcAnalysis(g_del)

result.valid % Indicates if the system robustly stable

result.performance % Bound on the worst case l2-induced norm for any admissible delta

Analyzing an uncertain system: Part II

But what if our uncertainty were time-varying? In other words, consider the (now time-varying) system:

In this case, we now make a time-varying, rather than time-invariant, delta:

% A Static, Linear, Time-Varying (SLTV) uncertainty

del_ltv = DeltaSltv('del', dims, lower_bound, upper_bound)

a_del_ltv = [0, 3/4 + del_ltv; 3/4 - del_ltv, 0];

g_del_ltv = toLft(a_del_ltv, b, c, d, timestep)

result = iqcAnalysis(g_del_ltv)

result.valid % This will be false, because the LFT is not robustly stable

result.performance % As this is not robustly stable, there is no upper-bound on the performance

An inability to conclude robust stability is expected. Consider the following admissible sample of del_ltv:

The resulting system is unstable. Look, for instance, at the system output given a non-zero initial condition and a zero input signal:

horizon_period = [0, 2]; % for a a 2-periodic system (0 initial horizon)

% Create a sample admissible time-varying parameter

delta_sample = toLft({-1/2 * eye(2), 1/2 * eye(2)}, horizon_period)

% Match horizon period of uncertain system to that of sampled delta

g_del_ltv_2periodic = g_del_ltv.matchHorizonPeriod(horizon_period);

% Insert sampled delta in place of the uncertainty 'del'

g_del_ltv_sampled = g_del_ltv_2periodic.sampleDeltas({'del'}, {delta_sample})

x0 = [1; 0]; % Initial condition of system

N = 31; % timesteps

t = [0:N - 1]; % time (only defined for clarity)

u = zeros(1, N); % zero input signal

y = simulate(g_del_ltv_sampled, u, t, x0);

% Simulate sampled system with given input and initial condition

figure

plot(y)

Here we see that the null result from IQC analysis on G_del_ltv was consistent with the system NOT being robustly stable. However, when we narrow the uncertainty set to only consist of time invariant uncertainties (holding all other properties the same), IQCs can leverage that additional information and correctly conclude robust stability.

Analyzing an uncertain system: Part III

A key advantage of IQC analysis is the broad library of IQCs that can characterize sets of uncertainties. In the last example, we saw that the uncertain system G_del was robustly stable if the uncertainty were time invariant, but NOT robustly stable if the uncertainty were time varying. But what if del were time-varying and it had a bound on the rate-of-change? There are IQCs to express that uncertainty set too.

upper_rate = 1/4;

lower_rate = -1/4;

% A Static, Linear, Time Varying, Rate Bounded uncertainty, with bounds and rates as given

del_ltv_rate = DeltaSltvRateBnd('del', dims, lower_bound, upper_bound, lower_rate, upper_rate);

a_del_ltv_rate = [0, 3/4 + del_ltv_rate; 3/4 - del_ltv_rate, 0];

g_del_ltv_rate = toLft(a_del_ltv_rate, b, c, d, timestep)

result = iqcAnalysis(g_del_ltv_rate)

result.valid

result.performance

Here we see that the resulting system is robustly stable. If we set the rate bounds to be zero, we recover the results pertaining to time-invariant uncertainties.

del_ltv_rate_zero = DeltaSltvRateBnd('del', dims, lower_bound, upper_bound, 0, 0);

a_del_ltv_rate_zero = [0, 3/4 + del_ltv_rate_zero; 3/4 - del_ltv_rate_zero, 0];

g_del_ltv_rate_zero = toLft(a_del_ltv_rate_zero, b, c, d, timestep)

result = iqcAnalysis(g_del_ltv_rate_zero)

result.valid

result.performance

Additionally, if we set the rate bounds to be too big, we recover the results pertaining to time-varying, arbitrarily fast uncertainties.

del_ltv_rate_fast = DeltaSltvRateBnd('del', dims, lower_bound, upper_bound, -1, 1);

a_del_ltv_rate_fast = [0, 3/4 + del_ltv_rate_fast; 3/4 - del_ltv_rate_fast, 0];

g_del_ltv_rate_fast = toLft(a_del_ltv_rate_fast, b, c, d, timestep)

result = iqcAnalysis(g_del_ltv_rate_fast)

result.valid

result.performance

The ability to flexibly and modularly express a variety of uncertainties is a core utility of IQC analysis. While the available IQCs in the literature are vast [1], [2], not all are included in iqcToolbox. These can be added to iqcToolbox as subclasses of the abstract Delta class (more on this in a separate section).

Reachability Analysis

The robust L2-induced performance metric asserts a bound on the ratio between the energies of a system's output signal and the system's input signal. This is a useful metric, but more intuitive ones may be asserted with IQCs. One such metric is the robust L2-to-Euclidean performance, which asserts a bound on the ratio between the Euclidean norm of a system's output (at a specific time instant) and the system's input signal energy.

The latter metric is capable of asserting a bound on a system's output at a specific time instant, given a bound on the input signal's energy. This can be measured with IQC analysis by creating and analyzing a finite-horizon system whose output is zero at any time except for the desired time instant.

Let's try this on one of the previous LFTs

final_time = 30; % Unit-less timesteps

g_del_ltv_reach = generateReachabilityLft(g_del_ltv, final_time)

Observe that the g_del_ltv_reach.horizon_period is now [final_time + 1, 1], instead of the former [0, 1]. This is because g_del_ltv_reach represents a finite-horizon system, which has zero output at k < final_time, normal output at k == final_time, and then has all zero matrices afterwards (in the 1-periodic final portion). If horizon_period is a novel concept to you, consider reading more about the Ulft class here. Let's analyze the system.

result = iqcAnalysis(g_del_ltv_reach)

result.valid

result.performance

Notice how the system is now robustly stable, even though g_del_ltv was not. This is because the final sequence of A-matrices for g_del_ltv are all zeros; since we only care about the system output at a specific time, the system dynamics don't matter after that final_time, and so are all set to zero for the remainder of time, which constitutes a stable system. This allows us to apply reachability analysis with IQCs on systems that are not stable.

These results now means that, if the L2-norm of the input signal equals 1, then y(final_time + 1) < result.performance. You can verify this yourself by generating input signals and simulating the system.

Incorporating Disturbances

IQCs can also be used to characterize sets of disturbances. Whereas the previous results assumed admissible disturbances could be any L2 signal, we can restrict the admissible disturbances to be a subset of L2. This changes our previous definition of performance to treat subsets of L2 signals:

For example, let's make a disturbance set which is non-zero ONLY after 15 timesteps.

horizon_period = g_del_ltv_reach.horizon_period;

dis = DisturbanceTimeWindow('after15', {[1]}, 15:(sum(horizon_period) - 1), horizon_period)

g_del_ltv_reach_dis = addDisturbance(g_del_ltv_reach, {dis})

NOTE:

The system time is zero-indexed. Hence, for a system whose horizon_period is [h, p], the maximum time instant to be specified can be h + p - 1. Also, any time-instants specified as non-zero AFTER h are repeated at each period. In this example, 17 time-instants are listed as non-zero (15th to 31st), but that number just pertains to the horizon and the FIRST period. In reality, an infinite number of time-instants are non-zero (15th and greater).

Now let's analyze this system whose disturbance is thus constrained.

result = iqcAnalysis(g_del_ltv_reach_dis)

result.performance

For this specific problem, the performance for the previous system would be equivalent to the performance of a system whose input begins at the 0th (rather than 15th) time-instant, and whose output is inspected at the 15th (rather than 30th) time-instant. Let's check that.

g_del_ltv_reach_15 = generateReachabilityLft(g_del_ltv, 15)

result = iqcAnalysis(g_del_ltv_reach_15)

result.performance

Incorporating uncertain initial conditions

So far we've discussed how model uncertaintites and disturbances can be incorporated in IQC analysis. We can also incorporate uncertain initial conditions and measure their impact on system performance. Consider our previous system, which may also have an initial condition within the unit sphere. Let's analyze the impact of that uncertain initial condition on the system output.

e_mat = eye(2); % Ellipse specified as {x in R^n | x' * e_mat * x <= 1}

uncertain_states = [true, true]; % Indicates first and second states are considered uncertain

options = AnalysisOptions('init_cond_ellipse', e_mat, 'init_cond_states', uncertain_states);

result = iqcAnalysis(g_del_ltv_reach, 'analysis_options', options)

Here we see new fields in the result variable that were previously NaN.

result.ellipse

result.state_amplification

Let's define what these results mean.

In our previous analyses, we assumed the initial conditions were zero, and inspected how disturbances could affect the system output. Now, we want to understand how the disturbances AND the initial conditions affect the system output.

Note:

Additional details on this updated metric is provided in [3].

Whereas the former performance metric bounded the system output energy only in terms of the system input energy, the updated performance metric bounds the system output energy with a combination of the system input energy, and the "state amplication" (lambda in the previous definition). In other words, result.state_amplification additionally contributes to the bound on the system output, and IQC analysis is now optimizing two decision variables, gamma (or result.performance), and lambda (or result.state_amplification).

To scalarize this optimization, you can set a weighting factor on the state_amplification term.

weight = 100;

options = AnalysisOptions('init_cond_ellipse', e_mat,...

'init_cond_states', uncertain_states,...

'scale_state_obj', weight);

A greater scale_state_obj places more emphasis on minimizing state_amplification. Here's what the paretol optimal front looks like for the previous system.

% Set a list of weights to solve for pareto optimal front

options.yalmip_settings.verbose = false;

weights = [0, logspace(-3, 6, 6)];

gamma = zeros(length(weights), 1);

lambda = gamma;

for i = 1:length(weights)

options.scale_state_obj = weights(i);

result = iqcAnalysis(g_del_ltv_reach, 'analysis_options', options);

gamma(i) = result.performance;

lambda(i) = result.state_amplification;

end

% Fix smallest gamma value, then find the smallest lambda value (find leftmost point on curve)

[dim_out, dim_in] = size(g_del_ltv_reach);

gam = gamma(1) * 1.01;

perf_smallest = PerformanceL2Induced('smallest',{1}, {1}, gam, g_del_ltv_reach.horizon_period);

g_del_ltv_reach_perf = addPerformance(g_del_ltv_reach, {perf_smallest});

options.scale_state_obj = 1e-4;

options.lmi_shift = 1e-7;

result = iqcAnalysis(g_del_ltv_reach_perf, 'analysis_options', options);

lambda(1) = result.state_amplification;

% Prioritize minimizing lambda value, with no cost to gamma value (find rightmost point on curve)

mult_perf = MultiplierL2Induced(g_del_ltv_reach.performance.performances{1}, dim_out, dim_in,...

'objective_scaling', 0);

% Last command created a specific multiplier for PerformanceL2Induced

result = iqcAnalysis(g_del_ltv_reach,...

'analysis_options', options,...

'multipliers_performance', mult_perf);

gamma(end + 1) = double(result.multipliers_performance.gain);

lambda(end + 1) = result.state_amplification;

% Plot values

figure

plot(gamma, lambda, '-*')

title('Robust $(E(\Xi),\mathcal{D})$-to-$\ell_2$-induced norm, $(\lambda, \gamma)$-values',...

'interpreter', 'latex')

xlabel('$\gamma$', 'Interpreter', 'latex')

ylabel('$\lambda$', 'Interpreter', 'latex')

If we know beforehand the bound on the input signal's L2-norm, we can use that to find the tightest bound on the output energy while incorporating the effect of uncertain initial conditions.

bound_l2_norm_in = 1;

mult_perf = MultiplierL2Induced(g_del_ltv_reach.performance.performances{1},...

dim_out, dim_in,...

'objective_scaling', bound_l2_norm_in^2);

options.scale_state_obj = 1;

result = iqcAnalysis(g_del_ltv_reach,...

'analysis_options', options,...

'multipliers_performance', mult_perf);

result.performance % Optimal gamma value (for given bound on input signal energy)

result.state_amplification % Optimal lambda values (for given bound on input signal energy)

opt_bound = sqrt(bound_l2_norm_in^2 * result.performance^2 + result.state_amplification^2)

Recalling the section on reachability analysis, this bound signifies that y(30) <= opt_bound.

Combining the disturbance specification capabilities from the previous section, let's find the tightest bound on the system output when there are no disturbances. Here, we'll set the disturbance to be non-zero after the finite horizon.

dis = DisturbanceTimeWindow('after30', {1}, 31, g_del_ltv_reach.horizon_period);

g_del_ltv_reach_30 = addDisturbance(g_del_ltv_reach, {dis})

perf_no_gain = PerformanceL2Induced('no_gain',{1}, {1}, 1e-4, g_del_ltv_reach.horizon_period);

g_del_ltv_reach_30 = addPerformance(g_del_ltv_reach_30, {perf_no_gain})

result = iqcAnalysis(g_del_ltv_reach_30, 'analysis_options', options)

result.performance

result.state_amplification

We see that this result is tight. Simulating the system with the following samplings:

yields an output at y(31) equal to result.state_amplification

x0

u = zeros(1, 31); % zero input

t = 0:30; % time starting from k = 0

y = simulate(g_del_ltv_sampled, u, t, x0);

figure

plot(y)

y(31)

Conclusion

Congratulations, you reached the end of this tutorial! By now you should have sufficient familiarity with the tools to generate LFTs, specify their uncertainties, specify disturbances, conduct reachability analysis, and incorporate uncertain initial conditions in your analysis. More in-depth tutorials are provided for generating LFTs, specifying uncertainties, disturbances, and performances, and conducting IQC analysis.

References

[1] Megretski, Alexandre, and Anders Rantzer. "System analysis via integral quadratic constraints." IEEE Transactions on Automatic Control 42, no. 6 (1997): 819-830.

[2] Veenman, Joost, Carsten W. Scherer, and Hakan Köroğlu. "Robust stability and performance analysis based on integral quadratic constraints." European Journal of Control 31 (2016): 1-32.

[3] Fry, J. Micah, Dany Abou Jaoude, and Mazen Farhood. "Robustness analysis of uncertain time‐varying systems using integral quadratic constraints with time‐varying multipliers." International Journal of Robust and Nonlinear Control 31, no. 3 (2021): 733-758.