MATLAB Examples

Simulink example with demo 'ab' in DST mode

This demonstration uses the 'ab' example of AUTO to demonstrate the use of Simulink with the Dynamical Systems Toolbox in the DST mode.

Contents

Open simulink model

The equations are defined in a Simulink model "ab.mdl".

open_system('ab');
snapnow;

Find the equilibrium of the system

The continuation run has to start from an equilibrium state, or a stable limit cycle condition. Most applications will start from an equilibrium condition. This equilibrium is also known as a quasi-steady state, or a trim condition. There are two ways to obtain such an equilibrium:

  1. Run the simulation for a sufficient period of time, until the derivatives (of the states of interest) are equal to zero. This is the easiest option.
  2. Use a trim routine to obtain a trim point.
clear all

Choose/guess initial conditions

PAR0=[0,14,2;0,14,2];
ts=[0,100];
U0=[0,0];
UT=[ts',PAR0];

options=simset('InitialState',U0);
[T,X,Y]=sim('ab',ts,options,UT);

Plot the results out. We can see that the state derivatives are zero. This is a trivial example seeing that we could have checked this by hand, but the method is applicable for more complex models.

h=plot(T,Y);
set(h,'LineWidth',2);
title('Time history for Simulink Demo: ab');
xlabel('Time (s)');
ylabel('U(1),U(2)');
snapnow;
close(gcf);

Adjust equations file to contain correct initial conditions We can now use the final states of the simulation as starting conditions for the continuation runs. If the derivatives for these states are not zero, the continuation runs will fail.

Create an auto object

a{1}=auto;

Define the function file and intial conditions

a{1}.s.FuncFileName='absimfunc';
a{1}.s.SimulinkModel='ab';
a{1}.s.Par0=PAR0(end,:);
a{1}.s.U0=X(end,:);
a{1}.s.Out0=Y(end,:);

Print function file to screen.

type(a{1}.s.FuncFileName);
function [f,o,dfdu,dfdp]=absimfunc(par,u,ijac)
%
% function file for demo ab
%
f=[0,0];
o=[];
dfdu=[];
dfdp=[];

% Call outputs first, otherwise derivatives calulated incorrectly. Not
% sure why.
o=ab(0,u,par(1:3),'outputs');
f=ab(0,u,par(1:3),'derivs');


Set the constants and run the stationary solutions

We can set the constants from the command line, but for ease of use we can convert the constants file from the original fortran code with the convertchc function, and then we can read the constants in from the m-file.

a{1}.c=cab1(a{1}.c);

Show format of set commands

type('cab1.m');
function c=cab1(c)
% cab1 - Constants file converted with function convertchc from c. format
%
%  Created by : ecoetzee
%  Date       : 15-Oct-2010 13:31:49
%

set(c,'Ndim',2,'Noutx',200,'Ips',1,'Irs',0,'Ilp',1);
set(c,'Icp',[1]);
set(c,'Ntst',50,'Ncol',4,'Iad',3,'Isp',1,'Isw',1,'Iplt',0,'Nbc',0,'Nint',0);
set(c,'Nmx',100,'Rl0',0,'Rl1',0.15,'A0',0,'A1',100);
set(c,'Npr',100,'Mxbf',10,'Iid',2,'Itmx',8,'Itnw',5,'Nwtn',3,'Jac',0);
set(c,'Epsl',1e-006,'Epsu',1e-006,'Epss',0.0001);
set(c,'Ds',0.01,'Dsmin',0.005,'Dsmax',0.05,'Iads',1);
set(c,'Thl',[10,0]);
set(c,'Thu',[]);
set(c,'Uzr',[]);

Print constants to screen.

a{1}.c
ans = 

  autoconstants handle

  Properties:
     Ndim: 2
    Noutx: 200
      Ips: 1
      Irs: 0
      Ilp: 1
      Icp: 1
     Ntst: 50
     Ncol: 4
      Iad: 3
      Isp: 1
      Isw: 1
     Iplt: 0
      Nbc: 0
     Nint: 0
      Nmx: 100
      Rl0: 0
      Rl1: 0.1500
       A0: 0
       A1: 100
      Npr: 100
     Mxbf: 10
      Iid: 2
     Itmx: 8
     Itnw: 5
     Nwtn: 3
      Jac: 0
     Epsl: 1.0000e-006
     Epsu: 1.0000e-006
     Epss: 1.0000e-004
       Ds: 0.0100
    Dsmin: 0.0050
    Dsmax: 0.0500
     Iads: 1
      Thl: [10 0]
      Thu: []
      Uzr: []


Run an equilibrium solution

Compute stationary solutions by using the runauto method in the auto object. The Simulink model will be opened if it is not open, and then compiled automatically before the analysis starts.

a{1}=runauto(a{1});
 
    --------------- DYNAMICAL SYSTEMS TOOLBOX ---------------------     
 
USER NAME      : ECOETZEE
DATE           : 26/10/2010 10:03:50
 
SIMULINK MODEL : ab
 
<
  BR    PT  TY  LAB      PAR(01)      L2-NORM         U(01)         U(02)
   1     1  EP    1   0.00000E+00   0.00000E+00   0.00000E+00   0.00000E+00
   1    33  LP    2   1.05739E-01   1.48439E+00   3.11023E-01   1.45144E+00
   1    70  LP    3   8.89318E-02   3.28824E+00   6.88982E-01   3.21525E+00
   1    90  HB    4   1.30900E-01   4.27187E+00   8.95080E-01   4.17704E+00
   1    92  EP    5   1.51242E-01   4.36975E+00   9.15589E-01   4.27275E+00

 Total Time    0.531E+00
>

Check contents written to objects

In the Fortran version of AUTO the results are written to the fort.7, fort.8 and fort.9 files. This is still an option in the DST mode, but it is not the default.

Print f7 contents to screen.

a{1}.f7
ans = 

  autof7

  Properties:
       Ibr: [92x1 double]
      Mtot: [92x1 double]
       Itp: [92x1 double]
       Lab: [92x1 double]
       Par: [92x1 double]
    L2norm: [92x1 double]
         U: [92x2 double]
       Out: [92x2 double]


Print f8 contents to screen.

a{1}.f8
ans = 

  autof8

  Properties:
       Ibr: [5x1 double]
      Mtot: [5x1 double]
       Itp: [5x1 double]
       Lab: [5x1 double]
      Nfpr: [5x1 double]
       Isw: [5x1 double]
      Ntpl: [5x1 double]
       Nar: [5x1 double]
    Nrowpr: [5x1 double]
      Ntst: [5x1 double]
      Ncol: [5x1 double]
     Nparx: [5x1 double]
      Ifpr: []
         T: [5x1 double]
        Tm: []
       Par: [5x36 double]
     Rldot: []
         U: [5x2 double]
       Ups: []
    Udotps: []


Plot stationary solutions

We have seen from the output to the MATLAB command window that there are two limit points and one Hopf bifurcation. We use the plautobj object to plot the results. Note the solid lines are used for stable equilibrium states, and broken lines for unstable equilibrium states. The limit points (LP) indicate a qualitative change in the stability, while the Hopf bifurcation (HB) indicates a transition from a steady state to a limit cycle.

p=plautobj;
set(p,'xLab','Par','yLab','L2-norm','axLim',[0,0.16,0,10]);
ploteq(p,a);
snapnow;

Compute periodic solutions

We want to do the continuation from the Hopf Bifurcation to determine the amplitude and frequency of the limit cycle as we change the parameter.

Copy the information in the first object to the second one. Set constants by cloning constants object. Remember that we are using handle classes and that we do not want them to act like pointers, hence use copy method from the autoconstants object.

a{2}=a{1};
a{2}.c=copy(a{2}.c);
a{2}.c=cab2(a{2}.c);

% Run continuation
a{2}=runauto(a{2});
 
    --------------- DYNAMICAL SYSTEMS TOOLBOX ---------------------     
 
USER NAME      : ECOETZEE
DATE           : 26/10/2010 10:03:52
 
SIMULINK MODEL : ab
 
<
  BR    PT  TY  LAB      PAR(01)      L2-NORM     MAX U(01)     MAX U(02)      PERIOD
   4    30        6   1.21347E-01   4.08243E+00   9.83711E-01   6.30390E+00   2.29801E+00
   4    60        7   1.18591E-01   3.76857E+00   9.97665E-01   8.25264E+00   3.66907E+00
   4    90        8   1.14794E-01   3.11323E+00   9.99585E-01   9.97959E+00   6.25524E+00
   4   120        9   1.06928E-01   2.51769E+00   9.99410E-01   9.62126E+00   9.53437E+00
   4   150       10   1.05558E-01   2.14343E+00   9.99125E-01   9.33018E+00   1.60574E+01
   4   180       11   1.05507E-01   1.91447E+00   9.99204E-01   9.29964E+00   2.90787E+01
   4   210       12   1.05507E-01   1.81310E+00   9.99102E-01   9.29900E+00   4.40753E+01
   4   240       13   1.05507E-01   1.76109E+00   9.99096E-01   9.29896E+00   5.90741E+01
   4   270       14   1.05507E-01   1.72940E+00   9.99089E-01   9.29853E+00   7.40735E+01
   4   300  EP   15   1.05507E-01   1.70806E+00   9.99124E-01   9.29874E+00   8.90733E+01

 Total Time    0.294E+03
>

You can see that the runs are a lot longer than the case where we used a Matlab function. This is because you make twice as many function evaluations due to a call to get the outputs, and then a call to get the derivatives. The output call needs to be done first otherwise the derivatives might be incorrect. Not sure why this is the case, but we have seen instances where this occurs.

Plot periodic solutions

Add the Maximum values of the Hopf-bifurcation.

ploteq(p,a);
snapnow;

Run simulation with slow ramp input

A continuation analysis is similar to running a simulation with a very slow ramp input. We will thus rerun the simulation and then overlay the results on top of the bifurcation diagram. In our previous example we plotted the L2-NORM. We will now work it out from time history data. The stable solution falls close to or on top of the blue branches, and then the oscillations appear in the area of point 4, where the Hopf-bifurcation occurs, indicating the onset of limit cycles.

Define ramp input

PAR0=[0,14,2;0.2,14,2];
ts=[0,20000];
U0=[0,0];
UT=[ts',PAR0];

Run simulation.

options=simset('InitialState',U0);
[T,X,Y]=sim('ab',ts,options,UT);

Overlay simulation onto bifurcation diagram

Plot the simulation on top of the bifurcation diagram.

PAR=interp1(ts,PAR0(:,1),T);
L2norm=sqrt(Y(:,1).^2+Y(:,2).^2);
plot(PAR,L2norm,'m');
snapnow;
close(gcf)

Two-parameter continuation, follow locus of lower limit point (Forward)

Trace out the locus of the lower limit point while stepping in a forward direction. We vary the first and the thrid parameters for the continuation, so restart from label 2 in this first continuation, hence we need to copy the restart file from the first run to the appropriate name. Create new object and set constants

a{3}=a{1};
a{3}.c=copy(a{3}.c);
a{3}.c=cab3(a{3}.c);

% Run continuation
a{3}=runauto(a{3});
 
    --------------- DYNAMICAL SYSTEMS TOOLBOX ---------------------     
 
USER NAME      : ECOETZEE
DATE           : 26/10/2010 10:08:54
 
SIMULINK MODEL : ab
 
<
  BR    PT  TY  LAB      PAR(01)      L2-NORM         U(01)         U(02)      PAR(03)
   2    27  LP    6   1.35335E-01   2.05992E+00   4.99604E-01   1.99842E+00   2.50000E+00
   2   100  EP    7   2.19708E-07   1.81989E+01   9.44978E-01   1.81744E+01  -2.72070E-01

 Total Time    0.616E+01
>

Two-parameter continuation, follow locus of lower limit point (Backward)

Trace out the locus of the lower limit point while stepping in a backward direction. We vary the first and the thrid parameters for the continuation, so restart from label 2 in this first continuation, hence we need to copy the restart file from the first run to the appropriate name. Create new object and set constants

a{4}=a{1};
a{4}.c=copy(a{4}.c);
a{4}.c=cab4(a{4}.c);

% Run continuation
a{4}=runauto(a{4});
 
    --------------- DYNAMICAL SYSTEMS TOOLBOX ---------------------     
 
USER NAME      : ECOETZEE
DATE           : 26/10/2010 10:09:00
 
SIMULINK MODEL : ab
 
<
  BR    PT  TY  LAB      PAR(01)      L2-NORM         U(01)         U(02)      PAR(03)
   2    36  EP    6  -3.51788E-03   9.90574E-01  -9.56303E-03   9.90528E-01  -1.13516E+00

 Total Time    0.219E+01
>

Two-parameter continuation, follow locus of hopf point (Backward)

Trace out the locus of the Hopf Bifurcation while stepping in a backward direction. We vary the first and the thrid parameters for the continuation, so restart from label 4 in this first continuation, hence we need to copy the restart file from the first run to the appropriate name. Create new object and set constants

a{5}=a{1};
a{5}.c=copy(a{5}.c);
a{5}.c=cab5(a{5}.c);

% Run continuation
a{5}=runauto(a{5});
 
    --------------- DYNAMICAL SYSTEMS TOOLBOX ---------------------     
 
USER NAME      : ECOETZEE
DATE           : 26/10/2010 10:09:03
 
SIMULINK MODEL : ab
 
<
  BR    PT  TY  LAB      PAR(01)      L2-NORM         U(01)         U(02)      PAR(03)
   4   100  EP    6   8.81006E-05   1.17440E+01   9.14609E-01   1.17083E+01   9.36293E-02

 Total Time    0.833E+01
>

Plot the locii of the limit points and Hopf bifurcation

Add the locii of the limit points and Hopf bifurcation

ploteq(p,a);
snapnow;
close(gcf);