## Contents

- Define the operator (T x)(s) = int_{a}^{s} x(r) dr on [a,b]=[0,1]
- Solve the LPI Top'*Top - gam<=0

    % DEMO2_volterra_operator_norm.m
    % See Chapter 11.2 of the manual for a description.
    %
    % This document illustrates how the norm of the Volterra integral operator
    % can be computed in PIETOOLS.

    % This example is also included in the paper (page 5, Demonstration 2)
    % link: https://arxiv.org/pdf/1910.01338.pdf

    % Volterra integral operator
    % T x (s) = int(x(\theta),d\theta,a,s)
    % Domain:  s,theta \in [a,b] = [0,1]

    % Optimization Problem (LPI)
    % min γ, such that
    % T^*T ≤ γ

    clc; clear;
    echo on
    %%%%%%%%%%%%%%%%%% Start Code Snippet %%%%%%%%%%%%%%%%%%

## Define the operator (T x)(s) = int_{a}^{s} x(r) dr on [a,b]=[0,1]

    a=0;    b=1;
    opvar Top;
    Top.R.R1 = 1;   Top.I = [a,b];

## Solve the LPI Top'*Top - gam<=0

First, define dpvar gam and set up an optimization problem

    vars = [Top.var1;Top.var2];     % Free vars in optimization problem (no optimization over these vars)
    dpvar gam;                      % Decision var in optimization problem (we will minimize gam)
    prob = sosprogram(vars,gam);

    % Next, set gam as objective function min{gam}
    prob = sossetobj(prob, gam);

    % Then, enforce the constraint Top'*Top-gam<=0
    opts.psatz = 1;                             % Allow Top'*Top-gam>0 outside of interval [a,b]
    prob = lpi_ineq(prob,-(Top'*Top-gam),opts); % lpi_ineq(prob,Q) enforces Q>=0

    % Finally, solve and retrieve the solution
    prob = sossolve(prob);
    operator_norm = sqrt(double(sosgetsol(prob,gam)));

    %%%%%%%%%%%%%%%%%% End Code Snippet %%%%%%%%%%%%%%%%%%
    echo off

    fprintf(['\n If successful, ',num2str(operator_norm),' is an upper bound on the norm of the Volterra integral operator.\n'])
    fprintf([' The exact operator norm of the Volterra integral operator is 2/pi=',num2str(2/pi),'.\n']);

    %%%%%%%%%%%%%%%%%% Start Code Snippet %%%%%%%%%%%%%%%%%%

    %%% Define the operator (T x)(s) = int_{a}^{s} x(r) dr on [a,b]=[0,1]
    a=0;    b=1;
    opvar Top;
    Top.R.R1 = 1;   Top.I = [a,b];

    %%% Solve the LPI Top'*Top - gam<=0
    % First, define dpvar gam and set up an optimization problem
    vars = [Top.var1;Top.var2];     % Free vars in optimization problem (no optimization over these vars)
    dpvar gam;                      % Decision var in optimization problem (we will minimize gam)
    prob = sosprogram(vars,gam);

    % Next, set gam as objective function min{gam}
    prob = sossetobj(prob, gam);

    % Then, enforce the constraint Top'*Top-gam<=0
    opts.psatz = 1;                             % Allow Top'*Top-gam>0 outside of interval [a,b]
    prob = lpi_ineq(prob,-(Top'*Top-gam),opts); % lpi_ineq(prob,Q) enforces Q>=0

    % Finally, solve and retrieve the solution
    prob = sossolve(prob);
    Size: 51  31

    The coefficient matrix is not full row rank, numerical problems may occur.
    SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
    Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
    Put 1 free variables in a quadratic cone
    eqs m = 31, order n = 13, dim = 53, blocks = 4
    nnz(A) = 130 + 0, nnz(ADA) = 933, nnz(L) = 482
     it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
      0 :            2.30E+00 0.000
      1 :   1.72E+00 5.88E-01 0.000 0.2562 0.9000 0.9000   0.49  1  1  2.4E+00
      2 :   7.92E-01 1.63E-01 0.000 0.2776 0.9000 0.9000   1.24  1  1  6.9E-01
      3 :   5.31E-01 3.84E-02 0.000 0.2349 0.9000 0.9000   1.43  1  1  1.4E-01
      4 :   4.71E-01 8.39E-03 0.000 0.2188 0.9000 0.9000   1.20  1  1  2.8E-02
      5 :   4.69E-01 2.16E-03 0.000 0.2568 0.9000 0.9000   0.83  1  1  7.3E-03
      6 :   4.69E-01 6.11E-04 0.000 0.2834 0.9000 0.9000   0.70  1  1  2.1E-03
      7 :   4.70E-01 1.73E-04 0.000 0.2831 0.9000 0.9000   0.73  1  1  7.0E-04
      8 :   4.70E-01 5.89E-05 0.000 0.3407 0.9000 0.9000   0.69  1  1  2.9E-04
      9 :   4.71E-01 2.55E-05 0.000 0.4320 0.9000 0.9000   0.62  1  1  1.6E-04
     10 :   4.71E-01 9.82E-06 0.000 0.3858 0.9000 0.9000   0.69  1  1  7.3E-05
     11 :   4.71E-01 4.71E-06 0.000 0.4795 0.9000 0.9000   0.52  1  1  4.8E-05
     12 :   4.71E-01 1.71E-06 0.000 0.3639 0.9000 0.9000   0.66  1  1  2.1E-05
     13 :   4.72E-01 7.35E-07 0.000 0.4292 0.9000 0.9000   0.49  1  1  1.2E-05
     14 :   4.72E-01 2.72E-07 0.000 0.3700 0.9000 0.9000   0.61  1  1  5.7E-06
     15 :   4.72E-01 1.18E-07 0.000 0.4338 0.9000 0.9000   0.44  1  1  3.6E-06
     16 :   4.72E-01 4.28E-08 0.000 0.3623 0.9000 0.9000   0.59  1  1  1.6E-06
     17 :   4.72E-01 1.78E-08 0.000 0.4167 0.9000 0.9000   0.43  1  1  9.8E-07
     18 :   4.72E-01 6.48E-09 0.000 0.3637 0.9000 0.9000   0.57  1  2  4.5E-07
     19 :   4.72E-01 2.71E-09 0.000 0.4187 0.9000 0.9000   0.42  1  2  2.8E-07
     20 :   4.72E-01 9.80E-10 0.000 0.3611 0.9000 0.9000   0.56  1  2  1.3E-07
     21 :   4.72E-01 1.63E-10 0.000 0.1663 0.7177 0.9000   0.41  1  2  7.7E-08
     22 :   4.72E-01 5.16E-11 0.000 0.3164 0.9000 0.9141   0.56  1  2  3.1E-08
     23 :   4.72E-01 1.90E-11 0.000 0.3694 0.9049 0.9000   0.50  2  2  1.6E-08
     24 :   4.72E-01 8.87E-12 0.000 0.4657 0.6531 0.9000   0.45  2  2  1.1E-08
     25 :   4.72E-01 3.06E-12 0.000 0.3450 0.9000 0.9197   0.62  2  2  4.7E-09

    iter seconds digits       c*x               b*y
     25      0.2   5.4  4.7194587777e-01  4.7194755947e-01
    |Ax-b| =   8.0e-09, [Ay-c]_+ =   3.0E-10, |x|=  6.2e+02, |y|=  4.3e+03

    Detailed timing (sec)
       Pre          IPM          Post
    1.800E-02    1.310E-01    0.000E+00
    Max-norms: ||b||=1, ||c|| = 1,
    Cholesky |add|=0, |skip| = 15, ||L.L|| = 1.11599e+08.

    Residual norm: 8.0042e-09

             iter: 25
        feasratio: 0.6151
             pinf: 0
             dinf: 0
           numerr: 0
           timing: [0.0180 0.1310 0]
          wallsec: 0.1490
           cpusec: 0.1875

    operator_norm = sqrt(double(sosgetsol(prob,gam)));

    %%%%%%%%%%%%%%%%%% End Code Snippet %%%%%%%%%%%%%%%%%%
    echo off

     If successful, 0.68698 is an upper bound on the norm of the Volterra integral operator.
     The exact operator norm of the Volterra integral operator is 2/pi=0.63662.