MATLAB: How to plot an envelope above and below the data plot using the Optimization Toolbox 6.1 (R2011b)

Optimization Toolbox

After I plot my data (y vs. x), how can I plot an envelope around it? I would like to plot one curve above and one below my data.

Best Answer

  • I have attached the MATLAB code file “plot_data_and_envelope.m” which uses the following function as a starting point for approximating your data:
    y = a*tanh(b*(x-c)) + d;
    The variable ‘p’ is a vector of the parameters [a b c d]. I used the following initial value for p:
    p_initial = [40 0.02 450 40];
    I then used the function FMINUNC in the Optimization Toolbox to carry out an unconstrained optimization to refine this initial value, and obtained the following answer:
    p_opt =
    41.9599 0.0153 430.6022 37.4995
    More information about FMINUNC is available by typing the following in the MATLAB Command Window:
    web([docroot,'/toolbox/optim/ug/fminunc.html'])
    You may access the same information on the following webpage:
    If you put a debugging breakpoint at line 15 of the code (either by using the icon for setting a breakpoint in the MATLAB code editor, or by using the “Debug” menu in the code editor), you can view the effects of the parameters by typing the following commands in the MATLAB command prompt:
    clf
    hold on
    plot(x, data)
    plot(x, f(p_initial, x), 'r')
    plot(x, f(p_opt, x), 'g')
    This shows that while the vector of initial values with ‘p_initial’ does produce a curve that matches the data fairly well, the curve produced using ‘p_opt’ produces a better match to the data. Therefore, I used the value of ‘p_opt’ as the new value for ‘p_initial’.
    The code then executes the two calls to FMINCON (as before, one of these calls includes the nonlinear inequality constraint for the upper bound, while the other uses the similar constraint for the lower bound).
    This example function ‘f’ and the corresponding initial values are meant to get you started; feel free to modify the curve and choice of initial values to suit your needs.
    I have attached a MATLAB code file “plot_envelope.m” which shows how to use the Optimization Toolbox and the FMINCON function, along with a nonlinear inequality constraint, to plot envelope curves. You can execute the code and view the plots. This example should get you started – feel free to modify the example for the equation you are plotting.
    The code first creates a sample curve. It is based on an exponential: 1 – exp(-t).
    A sinusoidal component is added to make the exponential oscillate, more so at the larger values of time. This is meant to roughly mimic the kind of data you are plotting. The code is as follows:
    t = linspace(0, 5, 1e3);
    data = 1 - exp(-t);
    data = data + 0.02*t.*sin(2*pi*3*t);
    The code works on the premise that the upper and lower envelopes can be obtained using a modified exponential, where the modification is in the form of a multiplicative parameter ‘c’, where ‘c’ would be less than 1 for the lower envelope, and greater than 1 for the upper envelope. The code for this function is as follows:
    function y = f(c,t)
    y = c*(1 - exp(-t));
    The idea is to put the data into an optimization solver FMINCON, which was chosen because it allows very flexible constraints to be imposed, which is what will be used to make the envelopes. When working with optimization solvers, it is necessary to define an objective function, which is here chosen as the sum of the squared error between the data and the fit-curve, as follows:
    obj_func = @(c, t, data) sum((data - f(c,t)).^2);
    The solvers also need an initial condition. I used a value of 2. You want to make this as close to the actual value as possible, to maximize the chance that the solver finds the correct answer.
    Next, the FMINCON solver is called to generate the upper bound, via the following code:
    c_upper = fmincon(@(c) obj_func(c,t,data), c_initial, ...
    [], [], [], [], [], [], ...
    @(c) nonlincon_upper_bound(c, t, data) )
    Notice that the nonlinear constraint function is also passed to the solver. This is defined as a separate function later in the code, as follows:
    function [v_ineq, v_eq] = nonlincon_upper_bound(c, t, data)
    v_ineq = max(data - f(c,t));
    v_eq = [];
    Similarly, another call to FMINCON includes the nonlinear constraint for finding the lower bound, which is also included as another function in the code file. It is very similar to the upper bound, except that the order of (data – f(c,t)) is swapped to be (f(c,t) – data) so that the envelope curve is below the data curve, rather than above it.
    The rest of the code plots the data curve, along with the upper and lower envelope curves. As I said, you will want to modify the function definition to suit your application.