#51 Frequency Response Plots using COLL/PO toolbox

[Jonathan Ehrmann]

Hello,

I'm trying to calculate frequency response plots of dynamical system (in my understanding continuation of periodic orbits). My first try was more or less successful using the COLL toolbox. Now, I'm wondering if the PO toolbox is the better choice. While studying the documentation the question came up if it's possible to use the excitation frequence instead of the excitation period as continuation parameter and release the same time 'po.period'.

[Harry Dankowicz]

Hello, Jonathan.

Thank you for your interest in using COCO.

The short answer to your question is as follows:

Preferred toolbox
For the continuation of periodic orbits in a smooth dynamical system, the 'po' toolbox is preferred to the 'coll' toolbox, since the former makes use of the latter and then appends additional constraints imposing periodicity and, as appropriate, non-degeneracy using a phase condition. The demo in po/examples/bistable shows how to use the 'po' toolbox for continuation of periodic orbits in a periodically forced Duffing oscillator.

Gluing conditions
In the periodically forced case, you must manually impose an additional constraint that establishes the appropriate relationship between the excitation period and the period of the periodic orbit. In the demo in po/examples/bistable, which continues periodic solutions of period equal to the excitation period, this is done using the command

prob = coco_add_glue(prob, 'glue', uidx(maps.T_idx), uidx(maps.p_idx(1)));

which imposes the constraint that the difference between the interval duration of the trajectory segment (i.e., the period) and the first problem parameter (here, the excitation period) equals 0.

An equivalent alternative is given by the command

 prob = coco_add_func(prob, 'glue', @(p,d,u) deal(d, u(1)-u(2)), [], 'zero', 'uidx', uidx([maps.T_idx, maps.p_idx(1)]));

If you define your vector field in terms of the excitation frequency omega, instead of the excitation period, then you would need to impose a constraint that the difference between the interval duration of the trajectory segment and 2*pi/omega equal 0. This can be accomplished with this command

prob = coco_add_func(prob, 'glue', @(p,d,u) deal(d, u(1)-2*pi/u(2)), [], 'zero', 'uidx', uidx([maps.T_idx, maps.p_idx(1)])); 

where it is assumed that the excitation frequency is the first problem parameter.

A longer answer will follow as a separate post.

I hope this helps.

/Harry

[Harry Dankowicz]

A longer answer is as follows:

Since COCO is an object-oriented development platform, the software is agnostic to the way you choose to construct your zero functions and monitor functions. The COCO toolboxes package a limited set of problem constructors for problems that are common in basic bifurcation analysis. These are each wrapped around one or several calls to the COCO core constructor coco_add_func.

As you work with particular classes of problems, you may realize that you repeat the same actions for many different problems of interest to you. To avoid such repetition, you may choose to encode such sequences of actions in a higher-level command. For example, instead of having to call ode_isol2po, coco_get_func_data, and coco_add_glue in sequence, as is done in the demo in the po/examples/bistable folder, you may define a new command that takes the appropriate arguments and executes this sequence.

The demo in the coll/examples/linode folder shows how to define a continuation problem using the 'coll' toolbox to perform continuation of periodic orbits of period 2*pifor a harmonically excited, damped linear oscillator under variations in the stiffness using either an autonomous or non-autonomous encoding of the vector field.

It is reasonably straightforward to edit this encoding to allow continuation under variations in the excitation period or excitation frequency. To do this, you would need to modify the definition of the vector field and the boundary conditions. It would also be instructive to implement the linear oscillator problem using the 'po' toolbox constructors. For an autonomous encoding, consider replacing harmonic terms with the time histories of the components of the Hopf normal form.

For frequency-response analysis of piecewise-smooth systems see, e.g., the po/examples/bangbang demo.

A few final words on terminology:
* The label 'po.period' is assigned by the 'po' toolbox constructors to a continuation parameter that tracks the interval duraction of the corresponding trajectory segment. The latter is a continuation variable that is part of the definition of the domain of the trajectory collocation zero problem. Continuation parameters are additional unknowns that track the values of monitor functions defined in terms of continuation variables. Continuation parameters may be fixed (inactive) or allowed to vary (active). For a non-autonomous encoding, 'po.period' is initially inactive and must be released in order to be allowed to vary. For an autonomous encoding, it is initially active and need to be exchanged with an inactive parameter in order to be fixed.
* In the po/examples/bistable demo, the excitation period appears as a problem parameter in the definition of the vector field. Problem parameters are continuation variables. The call to ode_isol2po additionally defines three continuation parameters with labels 'T', 'A', and 'd' that track the values of the problem parameters corresponding to the excitation period, excitation amplitude, and damping coefficient. These are initially inactive and, therefore, need to be released, as appropriate, in order to be allowed to vary. In the call to coco, 'po.period' and 'T' are included in order to allow both to vary, since the encoding is non-autonomous. Allowing only one to vary would be a violation of the assumption that the solution manifold is one-dimensional.

I hope that helps.

/Harry

[Kola]

Hello Harry, In order to better understand the functionality of the 'po' toolbox, I attempted to reproduce the results of the linode example by encoding the system as an autonomous set of equations. However, when I run the continuation problem with variations in the stiffness of the oscillator, I don't find a solution manifold as the COCO output under type shows MX.

One step I took towards troubleshooting this issue was verifying that the Initial trajectory from the autonomous system that I was feeding to the coco problem was the same as the initial trajectory generated for the non-autonomous system. For comparison, the time history of x1 and x2 for both the autonomous and nonautonomous systems are shown in the figures named Timeseries_x1 and Timeseries_x2 respectively.

I am not sure in which direction to look to trouble shoot the issue of not landing on a solution manifold. Do you have any recommendations on what I need to pay more attention to in order to set up the autonomous system so that COCO gives me a solution manifold similar to what was produced for the nonautonomous linode example ?

I have attached the MATLAB script and function files I used for the COCO call in the attached folder.

[Harry Dankowicz]

Thank you for your question.

The po toolbox constructors build problems associated with periodic trajectory segments, i.e., those for which the two end points coincide. The hspo toolbox constructors can handle other types of relationships between the end points of successive segments, including jumps.

In the non-autonomous encoding of the linear oscillator, the two state variables are identical at the two end points, and so po works fine. But if you augment the state to include a phase variable, then this is not a periodic quantity (indeed, it changes by a distance 2*pi between the end points. For this, you need to use hspo.

Alternatively, you can just use coll directly, as is shown in the linode demo in the coll toolbox folder. Please have a look and see if this helps.

Kind regards,

Harry

[Harry Dankowicz]

One further observation: however you choose to define your vector field, you can always perform the conversion between excitation frequency and excitation period in your post-continuation analysis, e.g., when graphing the frequency-response diagram using coco_plot_bd or your own plotting routines.

[Jonathan Ehrmann]

Hello Harry,

thanks for this extremly detailed and fast answer! I'll give my best to get along with my problem applying your tips.

Have a nice weekend, Jonathan

[Jonathan Ehrmann]

Hey Harry,

my further question is how to understand the part @(p,d,u) deal(d, u(1)-2*pi/u(2)). Which function is meant with @(p,d,u)? I get that we somehow have to define u(1)-2*pi/u(2)=0.

[Harry Dankowicz]

Hi Jonathan,

This is Matlab syntax for an anonymous function of three arguments (p, d, and u) that returns two arguments, one of which evaluates to the input argument d and one of which evaluates to the expression u(1)-2*pi/u(2). Look at the Matlab manual for more information about anonymous functions.

The reason to use this particular syntax here is the requirement that zero and monitor functions added using coco_add_func must be of this form. An alternative is to provide a function handle to a named function (i.e., one encoded in a variable or in an m-file). That function would also need to take three input arguments and return two (as described above). You would need to look in the core tutorial or Recipes for Continuation to learn more about this.

I hope that helps.

/Harry

[Jonathan Ehrmann]

Hi Harry,

it took a while for me to have time to go back to that topic... First, thanks for your last tip, I managed to implement those. Nevertheless, for some "more complicated" parameters in my system I get the following warning by Matlab during continuation:

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  7.252075e-17.

I'm not sure how to deal with this warning. I have the feeling that the results are also not correct. Do you have any tips for me?

Thanks in advance,
Jonathan

[Harry Dankowicz]

Hello, Jonathan.

This is difficult to debug in the absence of your actual code. My gut guess is that you have made errors in the encoding of derivatives. That's where I usually start, rather than to assume that the error is due to something inherently remarkable about the problem at hand.

In case you are able to make use of the symcoco routines for automatic generation of vector fields and their derivatives (see symcoco-doc.pdf), consider comparing the results.

Best,

Harry

[Jonathan Ehrmann]

Thanks a lot, I think there was a problem with the derivatives!

Maybe a general question regarding that topic: is it sufficient to just let Coco calculate the derivatives numerically (easiest way to implement in my eyes) or is it better to always give analytical derivatives if possible?

Best, Jonathan

[Harry Dankowicz]

Hi Jonathan,

Your question is both general to numerical methods and specific to COCO. For the general answer, I refer you to a multitude of other sources that explore the computational cost and benefits of explicit derivatives. For COCO, and given the use of derivatives in detecting special points (e.g., bifurcations), my general practice is to encode derivatives of vector fields, at least to first order and in many cases to higher order, including directional derivatives, as appropriate. The symcoco package does this for you, but also incurs some computational overhead that would not be present if you encoded the derivatives directly according to some optimized execution sequence.

If you intend to use the adjoint toolbox constructors, explicit first-order derivatives are required and explicit second-order derivatives are useful.

In all cases, vectorized implementations are absolutely preferable. Indeed, much of the effort in encoding toolbox constructors for particular classes of problems (and, in particular for constrained design optimization with COCO) is spent on encoding vectorized forms of the problem derivatives.

Now, during development, you may find it convenient to first explore a problem (when possible) with numerical derivatives. You will typically see a significant speed-up when you then encode explicit derivatives, but you will also have the advantage of knowing what to expect, albeit at some cost of execution time. This may avoid the problem of having warnings of unknown origin thrown during execution due to incorrect encodings.

Finally, for COCO-compatible problems constructed with coco_add_func it is possible to include a flag for checking the accuracy of the Jacobian encoding (see the help for coco_add_func for details). At present, there is no such flag for the coll, ep, or po constructors, although it would not be difficult to create one. You can also use symcoco to check your own encoding, but that seems a bit superfluous, given that you may as well just use symcoco from the get go if you have access to it (at least for most medium-sized problems).

Best,

Harry

[Jonathan Ehrmann]

Hi Harry,

that is helpful, thank you!

Best, Jonathan

[Ali Kogani]

Hi Harry

I am using PO toolbox to calculate the frequency response of a vibration system with nonlinear stiffness. I have one problem that I do not know how to solve.

In the first step, I calculate the frequency response (just like the example for Duffing in the tutorials). However, I do not know how to access the time response (orbits) in each step. I need them because I need to apply constraints on the amplitude and/or the phase of the time responses while relaxing another parameter in the system.

Best,
Ali

[Harry Dankowicz]

Dear Ali,

Thank you for reaching out and for contributing to this thread.

I am wondering if the information in the Getting StartedwithCOCO tutorial document and demos (new with the August release of COCO) may be helpful to you. For example, Section 7.3 shows how you can monitor properties of periodic orbits during continuation and Section 7.4 shows how to track and constraint orbit maxima. There are other ways to do this as well, but perhaps this material is a good start.

Best,

Harry

[Ali Kogani]

Dear Harry,

Thank you for your fast response

I looked at section 7.3 and 7.4 as you said. I saw that it uses "po_add_bddat". I tried it with my problem and I managed to monitor the orbits. Thank you for your guidance.

However, I have used coco for some time now and I have always used "coco_add_func" for all my monitor and zero (constraints) functions. Since I know how to apply constraints using "coco_add_func", I was wondering if there is a way to feed "xbp" into "coco_add_func" (type zero) and define my constraint functions in there and apply them.

Best,
Ali

[Harry Dankowicz]

The function I had in mind you using was po_add_func and I suspect that's what you, in fact, did. This function is a wrapper to coco_add_func that also includes a preceding call to coco_get_func_data. Together these commands extract required information from the associated 'coll' instance and then impose a further constraint on the corresponding continuation variables. The reason for introducing po_add_func as part of the GettingStartedwithCOCO tutorial was to avoid the need to require i) explicit usage of coco_get_func_data and ii) an understanding of the way in which the 'coll' and 'po' toolboxes store data.

However, since you are comfortable working with such core commands, I recommend that you have a look at the source code for po_add_func as this should guide you in implementing a zero function that accomplishes the outcome you seek, as a special case of the more general functionality provided by po_add_func.

Kind regards,

Harry

[Ali Kogani]

Dear Harry

I truly appreciate your responses.
I got it working!

Best,
Ali

[Bo]

Dear Harry,

I have a problem plotting the frequency amplitude curve. I am following the example for plotting the frequency amplitude curve of bistable, the following codes are constructed. Unfortunately, MX is displayed. I do not see in the documentation how to adjust the continuation settings to avoid MX.

The other equation is about the frequency amplitude curve in relation to the excitation frequency Omega https://sourceforge.net/p/cocotools/tickets/51/#021b. Are there details on how to plot this?

Best Regards,

Bo

[Harry Dankowicz]

Dear Bo,

Thank you for your question.

The MX label indicates a failure of the algorithm to converge to a solution within the desired tolerance. This is a warning to you, not necessarily that there is anything wrong with your continuation settings, but perhaps that there may be something amiss or at least needing attention with your problem.

For example, I note that the initial value for the first problem parameter is (line 4 of demo_po.m) 2*pi/7. From line 5 of gen_sym_po.m, I conclude that p1 is the period of the excitation. If you are looking for a period-1 periodic orbit, you would then generate an initial solution guess over a time interval 2*pi/7 in duration. Since you have chosen a duration of 2*pi it appears that you are looking for a period-7 solution. That seems a bit unusual to me, but perhaps that is what you intended.

Next, when I graph different components of your simulated solution, I notice a great deal of variation. For example, the 9th component of your state vector varies between -0.0000004 and 0.0000004, which is a very small range considering that other components of your state vector vary between -0.006 and 0.006, say. Such dramatic differences in size are challenging for numerical routines that rely on tolerances to determine convergence. This includes the numerical integration routine that you use to generate the initial solution guess.

More important, several of the state vector components appear to exhibit significant oscillatory behavior even during a single period of excitation. As a result, a fine mesh is required to capture the detailed dynamics.

So, here's what I just did to get your code running. I changed line 3 of demo_po.m to read

funcs = {F(''),F('x'),F('p'),F('t')};

I then changed line 6 to read

[t0,x0]=ode45(@(t,x)funcs{1}(t,x,p0),[0 2*pi/7],x0(end,:)'); 

Finally, I added the following after line 12:

prob = coco_set(prob, 'coll', 'NTST', 100);

I hope that helps.

As for your second question, it appears to me that the line

coco_plot_bd('freq_resp', 'T', @(T) 2*pi./T, '||po.orb.x||_{L_2[0,T]}')

in your code already plots the solution norm against the excitation frequency.

Harry

[Bo]

Dear Harry,

Thank you for your reply.

I have changed the codes based on your comments. I have a few questions I would like to ask for your help.

1. In the line 23 of demo_po.m,
   cont_args = {1, {'po.period' 'T'}, [2*pi/200 2*pi/2]};
   it appears that the period of the system is limited between 2 pi/200 and 2 pi/2 (gen_sym_po.m and demo_po.m) . The result is shown in Figure 1.
   If I change the following to read
   p0=[2*pi/20; 1];
[~,x0]=ode45(@(t,x)funcs{1}(t,x,p0),[0 20*pi],[0.01; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]);
[t0,x0]=ode45(@(t,x)funcs{1}(t,x,p0),[0 2*pi/20],x0(end,:)');
cont_args = {1, {'po.period' 'T'}, [2*pi/20 2*pi/2]};
   the curve is plotted from omega=20 (gen_sym_po_1.m and demo_po_1.m) . So if I want to plot the frequency amplitude curve of the system with the jumping behavior, I have to use gen_sym_po.m and gen_sym_po_1.m. In fact, the period of each point on the curve varies between 2 pi/20 and 2 pi/2. Why can I not plot the curve using only the code gen_sym_po.m?

2. As you said, I can plot the frequency amplitude curve using
   coco_plot_bd('freq_resp', 'T', @(T) 2*pi./T, '||po.orb.x||_{L_2[0,T]}')
   If I change the state vector of the system without the excitation period T instead of the excitation frequency omega (gen_sym_po_2.m), how should I change the code dem_po.m? I still do not know what the parameter u means in the following code.
   prob = coco_add_func(prob, 'glue', @(p,d,u) deal(d, u(1)-2*pi/u(2)), [], 'zero', 'uidx', uidx([maps.T_idx, maps.p_idx(1)]));

Best Regards,

Bo

[Harry Dankowicz]

Hi Bo,

As I read your questions, I think it is important to separate out different components of the code.

1. Let's start with

coco_plot_bd('freq_resp', 'T', @(T) 2*pi./T, '||po.orb.x||_{L_2[0,T]}')

This command is executed once a run of the coco entry-point function is complete. It assumes that the run was named 'freq_resp' and that there was a continuation parameter named 'T' (or that a user-defined function responding to the bddat signal created a column with this header). It also assumes that the continuation problem included a periodic orbit object with identifier 'po.orb' so that the output data array includes a column with header '||po.orb.x||_{L_2[0,T]}'. The command uses the data in the column with header 'T' to calculate a column of data given by 2*pi/'T' and assigns this to the x-coordinate of each point in the diagram. If you replaced this line with

coco_plot_bd('freq_resp', 'T', @(T) T, '||po.orb.x||_{L_2[0,T]}')

or simply

coco_plot_bd('freq_resp', 'T', '||po.orb.x||_{L_2[0,T]}')

the x-coordinate would just be the value of 'T'.

Whether any of this makes sense depends on the meaning of the data in the T column, which is defined by the user, not by COCO.

1. Next, let's look at the command

prob = coco_add_func(prob, 'glue', @(p,d,u) deal(d, u(1)-2*pi/u(2)), [], 'zero', 'uidx', uidx([maps.T_idx, maps.p_idx(1)]));

The core constructor coco_add_func is here used to add a zero function to the continuation problem (note the 'zero' flag in the command). The function identifier is chosen as 'glue'.

The function depends on the elements of the vector of previously defined continuation variables indexed by the integers uidx([maps.T_idx, maps.p_idx(1)]). Execute this expression on the command line to see what the integers happen to be. These continuation variables correspond to the duration of the trajectory segment (i.e., the period of the solution) and the first problem parameter (which is defined by your vector field).

The function has no function data, thus the empty bracket before 'zero'.

The expression @(p,d,u) deal(d, u(1)-2*pi/u(2)) is an anonymous Matlab function (you can read more about such functions in the Matlab manual) that depends on three arguments (p, d, and u) and that takes the argument u (which contains two components, namely the period of the solution and the first problem parameter) and returns d and u(1)-2*pi/u(2). Since 'glue' here identifies a zero function, it follows that

u(1)-2*pi/u(2)=0

Whether this makes sense or not depends on the meaning of the first problem parameter. If this parameter is the excitation frequency, then the zero function imposes the constraint that the period of the solution equal the period of excitation.

1. The range [2*pi/200 2*pi/2] is the computational domain. You are correct that this means that continuation will proceed only as long as the period of the periodic solution (not necessarily the excitation period, since the former could be a multiple of the latter) remains within these bounds. If you change the line to

cont_args = {1, {'po.period' 'T'}, {[],[2*pi/200 2*pi/2]}};

then the computational domain is defined by the value of the continuation parameter 'T' remaining with these bounds. If you are looking at a period-one solution, then the result will be the same. Since the excitation frequency is given by 2*pi/T, where T is the period of excitation, whose value is tracked by T, the excitation frequency will remain between 2 and 200.

The code in gen_sym_po.m has no effect on the result of your computation. This code is just used to generate encodings of the vector field and its derivatives. It knows nothing about the problem you are trying to solve. There should be no reason to change this code just because you plan to perform a different analysis. Decide how you want to encode your vector field and then stick to it.

If you want to change the range of excitation frequencies, you can just change the limits on the excitation period. Alternatively (as I think you are doing in your second question), you can introduce another continuation parameter (called omega) that tracks the value 2*pi/p(1). Then you can define the computational domain directly in terms of omega .

I hope that helps.

Harry

[Bo]

Dear Harry,

Thank you for your reply.

I am not sure I understand what you mean.

1. why can we not calculate the entire frequency response curve at conce?
   Reason 1: If I take the excitation period T=2pi/5 as the initial parameter, so that in the jump region of the frequency response curve, the actual excitation period T may be not continue with the frequency response curve (the excitation period has changed, not continued). Instead, we need to constrain the curve from the other point, such as T=2pi/20. Finally, we can plot the whole response curve in the region of [2, 20] for the excitation frequency omega.
   Reason 2: From the command window, we can see that the actual count is 100. Is that why we can not calculate the entire frequency response curve at conce? If so, where are the calculated points set? (in other softwares, an additional set of calculated points is required to analyze the frequency response curve).

2. I have changed the parameter T to the excitation frequency omega. But the 'MX' still appears in demo-po_1.

Best Regards,

Bo

[Harry Dankowicz]

Hi Bo,

Continuation proceeds as long as solutions remain inside the computational domain AND as long as the maximum number of points (PtMX) has not been reached. The default value for PtMX is 100 in both directions along a one-dimensional solution manifold. You can change this to other values, including different values in each of the two directions, by using a command like

prob = coco_set(prob, 'cont', 'PtMX', [200 50]);

As to how much of the solution manifold that is covered by a particular setting for PTMX will also depend on the size of individual steps between solutions. I encourage you to look at the log file associated with a particular run to see whether the step size is maxing out and, if so, whether you can increase the maximable allowable step size to provide for larger steps between solutions points.

I was able to get your original code to run without any instances of 'MX' but haven't tried your most recent code. I encourage you to assess whether your initial solution guess (orbit and parameters) is, in fact, a good approximation for the problem you are trying to solve.

Best,

Harry