Using the Worstcase Solver - Demo 1

The worstcase solver is used to find the induced L2-to-L2 gain of a four-state nonlinear system.

Timothy J. Wheeler
Dept. of Mechanical Engineering
University of California, Berkeley

Contents

Introduction.

Consider a dynamic system of the form

where x(0)=0. Given positive scalars B and T, the goal is to maximize

subject to the constraint

Note: we only consider inputs and outputs defined on the interval [0,T].

System parameters.

This system is parameterized by the following constants:

lam = 1;
PL = 1;
gammaX = 1;
gammaR = 1;
A = 0.8;
tau = 1;
K0x = (-1/tau - A)/lam;
K0r = (1/tau)/lam;

Create a model of the system.

First, polynomial variables are created using the pvar command. Then, these variables are used to define the functions f and g, which are also polynomial variables.

pvar x1 xm zx zr r w
states = [x1;xm;zx;zr];
inputs = [r;w];

f(1,1) = A*x1 + lam*((zx+K0x)*x1 + (zr+K0r)*r) + w;
f(2,1) = (1/tau)*(-xm+r);
f(3,1) = -gammaX*x1*(x1-xm)*PL;
f(4,1) = -gammaR*r*(x1-xm)*PL;

g = ((zx+K0x)*x1 + (zr+K0r)*r) + w;

Then, a polysys object is created from the polynomials f and g.

sys = polysys(f,g,states,inputs);

The polynomial objects states and inputs specify the ordering of the variables. That is, by setting states(1) = x1, we specify that f(1) is the time derivative of x1.

Optimization parameters.

Use the following values for the optimization parameters (defined above):

T = 10;
B = 3;

The time vector t specifies the time window (T = t(end)) and the points at which the system trajectory is computed.

t = linspace(0,T,100)';

Set options for worstcase solver.

Create a wcoptions object that contains the default options.

opt = wcoptions();

Specify the maximum number of iterations and which ODE solver to use.

opt.MaxIter = 50;
opt.ODESolver = 'ode45';

Tell the solver to display a text summary of each iteration.

opt.PlotProgress = 'none';

Specify the optimization objective, and the bound on the input.

opt.Objective = 'L2';
opt.InputL2Norm = B;

Find worstcase input.

[tOut,x,y,u,eNorm] = worstcase(sys,t,opt);

Simulate with worstcase input.

We can only compute the worstcase input over a finite interval of time [0,T]. However, any response of the system that occurs after the input is "shut off" (i.e., u(t) = 0 for t > T) should contribute to our objective. Hence, we compute a more accurate value of the objective by continuing the simulation from the end of the previous trajectory with no input:

[te,xe,ye] = sim(sys,tOut,x(end,:)');
td = [tOut;tOut(2:end)+max(tOut)];
yd = [y;ye(2:end)];

The objective value over [0,T] is

eNorm
eNorm =
    4.7436

The objective value over [0,2T] is

eNormd = get2norm(yd,td)
eNormd =
    4.9622

Display results.

fprintf( 'The L2-to-L2 gain is %f\n', eNormd/B );

figure;
plot(tOut,u)
xlabel('Time, t')
ylabel('Input, u(t)')
title('Worst case input.')

figure;
plot(td,yd)
xlabel('Time, t')
ylabel('Output, y(t)')
title('Worst case output over extended time interval.')
The L2-to-L2 gain is 1.654050

Specifying a starting point.

By default, the worstcase solver starts with a constant input and then searches for a better input. Since this problem is nonconvex, this search may get "stuck" at a local optimum. We can help the solver by specifying a sensible starting point that is known to exhibit a large output.

load demo1_badInput
u0 = B * ubad/get2norm(ubad,tbad);
opt.InitialInput = u0;

Run solver again.

[tOut,x,y,u,eNorm] = worstcase(sys,t,opt);

Extend this simulation.

[te,xe,ye] = sim(sys,tOut,x(end,:)');
td = [tOut;tOut(2:end)+max(tOut)];
yd = [y;ye(2:end)];

The objective value over [0,T] is

eNorm
eNorm =
    5.0020

The objective value over [0,2T] is

eNormd = get2norm(yd,td)
eNormd =
    5.0029

Note that we achieve a larger value of the objective when we start the solver at u0.

Display new results.

fprintf( 'The L2-to-L2 gain is %f\n', eNormd/B );

figure;
plot(tOut,u)
xlabel('Time, t')
ylabel('Input, u(t)')
title('Worst case input.')

figure;
plot(td,yd)
xlabel('Time, t')
ylabel('Output, y(t)')
title('Worst case output over extended time interval.')
The L2-to-L2 gain is 1.667635