Plot function within optimisation algorithm does not work with subplots -

Hello, I am trying to visualise optimisation progress, where the mathematical model is fitted onto measured data. By using only one set of axes, the function works with no problems. Whenever I try to implement subplots to plot multiple values into one figure, I get the following error:

Error using optim.problemdef.OptimizationProblem/solve - Invalid or deleted object.

I am using lsqnonlin function, the options of which are defined by the following:

lsqoptions = optimoptions("lsqnonlin",...
                          Display = "iter-detailed",...
                          MaxIterations = 1.5e4,...
                          MaxFunctionEvaluations = 1e5,...
                          FunctionTolerance = 1e-18,...
                          StepTolerance = 1e-18...
                          );

lsqoptions.PlotFcn = @(x, optimValues, flag) plotprogress(x, optimValues, flag, ipe, plotting, ord, freq, stor_mod, loss_mod, dyn_mod, E);

Within the plotprogress function,I define variables to be plotted along with several other quantities. When trying to plot with subplots, the plot function looks like this:

function out = plotprogress(varargin)
% extracting variables from varargin....

% plotting
ax = subplot(3,1,3);
plotting.count = 1;
for i = 0:18:73
    semilogx(ax,freq(i+1:i+18),dyn_mod(i+1:i+18),'.','LineWidth',2,'MarkerSize',25,'Color',plotting.dotcol(plotting.count,:)); hold on; grid on;
    plotting.count = plotting.count+1;
end
semilogx(ax,freq,dyn_num,'LineWidth',2,Color=ipe.darkred)
hold off;
legend('20C','10C','0C','-10C','-20C','fit','Location','northwest')
xlabel('f_{master}')
ylabel('Dynamic modulus master')

out = false;
end

Every variable is defined, so this should not be a problem. Whenever I change the ax = subplot(3,1,3); into ax = gca; the plot function seems to be working with no problems whatsoever.

[Jan]

Please post a copy of the complete message, which reveals, in which line the error occur. Sharing this important inforamtion with the readers ist a good idea.

[Ondrej Farkas]

I will share the code as a whole, I guess it will be easier.

clc
clear
close all

% ipe colors
ipe.navy = [0 0 0.502];
ipe.darkred = [0.545 0 0];
ipe.orange = [1 0.647 0];
ipe.seagreen = [0.18 0.545 0.341];
ipe.darkmagenta = [0.545 0 0.545];

load("dma_mastercurve.mat")
dma_data = table2array(tabledata);
freq = dma_data(:,10);
stor_mod = dma_data(:,12);
loss_mod = dma_data(:,13);
dyn_mod = dma_data(:,14);
tan_del = dma_data(:,15);

%% Definition of the optimisation process

% Equilibrium stress E [MPa], for MR2 E = 2*(C10 + C01)
E = 0.41888;

% Order of Prony series model
ord = 5;

% Optimisation variables
mu = optimvar('mu',ord,'LowerBound',0);
tau = optimvar('tau',ord,'LowerBound',0);

% Initial conditions
ini.mu = 0.5*ones(1,ord);

% Prony series in frequency domain
stor_num = E;
loss_num = 0;
for i=1:ord
    stor_num = stor_num + mu(i) * ( ((freq * tau(i)).^2) ./ (1 + ((freq * tau(i)).^2)) );
    loss_num = loss_num + mu(i) * ( (freq * tau(i)) ./ (1 + ((freq * tau(i)).^2)) );
    ini.tau(i) = 10^(i-1);
end

% Loss factor
tan_num = loss_num ./ stor_num;

% Objective function based on moduli
obj_fun = sum((stor_num - stor_mod).^2 + (loss_num - loss_mod).^2 + (tan_num - tan_del).^2);


% Optimisation problem
lsqproblem = optimproblem("Objective",obj_fun);

% Optimisation options
lsqoptions = optimoptions("lsqnonlin",...
                          Display = "iter-detailed",...
                          MaxIterations = 1.5e4,...
                          MaxFunctionEvaluations = 1e5,...
                          FunctionTolerance = 1e-18,...
                          StepTolerance = 1e-18...
                          );

lsqoptions.PlotFcn = @(x, optimValues, flag) plotprogress(x, optimValues, flag, ipe, plotting, ord, freq, stor_mod, loss_mod, dyn_mod, E);
% lsqoptions.PlotFcn = 'optimplotfval';

tic
[sol,fval] = solve(lsqproblem,ini,Options=lsqoptions);
toc

%% Plot Function
function out = plotprogress(varargin)

% Extracting function variables
E = varargin{end};
dyn_mod = varargin{end-1};
loss_mod = varargin{end-2};
stor_mod = varargin{end-3};
freq = varargin{end-4};
ord = varargin{end-5};
mu = varargin{1}(1:ord);
tau = varargin{1}(ord+1:end);
plotting = varargin{end-6};
ipe = varargin{end-7};

% Plotting partial solution
stor_num = E;
loss_num = 0;
for i=1:ord
    stor_num = stor_num + mu(i) * ( ((freq * tau(i)).^2) ./ (1 + ((freq * tau(i)).^2)) );
    loss_num = loss_num + mu(i) * ( (freq * tau(i)) ./ (1 + ((freq * tau(i)).^2)) );
end
dyn_num = sqrt(stor_num.^2 + loss_num.^2);

plotting.count = 1;
f = gcf;
set(f, 'Color', [1 1 1]);
% ax = gca;


% Storage modulus
subplot(3,1,1)
for i = 0:18:73
    semilogx(freq(i+1:i+18),stor_mod(i+1:i+18),'.','LineWidth',2,'MarkerSize',25,'Color',plotting.dotcol(plotting.count,:)); hold on; grid on;
    plotting.count = plotting.count+1;
end
semilogx(freq,stor_num,'LineWidth',2,Color=ipe.darkred)
hold off;
legend('20C','10C','0C','-10C','-20C','fit','Location','northwest')
xlabel('f_{master}')
ylabel('Storage modulus master')

% Loss modulus
subplot(3,1,2);
plotting.count = 1;
for i = 0:18:73
    semilogx(freq(i+1:i+18),loss_mod(i+1:i+18),'.','LineWidth',2,'MarkerSize',25,'Color',plotting.dotcol(plotting.count,:)); hold on; grid on;
    plotting.count = plotting.count+1;
end
semilogx(freq,loss_num,'LineWidth',2,Color=ipe.darkred)
hold off;
legend('20C','10C','0C','-10C','-20C','fit','Location','northwest')
xlabel('f_{master}')
ylabel('Loss modulus master')

% Dynamic modulus
subplot(3,1,3);
plotting.count = 1;
for i = 0:18:73
    semilogx(freq(i+1:i+18),dyn_mod(i+1:i+18),'.','LineWidth',2,'MarkerSize',25,'Color',plotting.dotcol(plotting.count,:)); hold on; grid on;
    plotting.count = plotting.count+1;
end
semilogx(freq,dyn_num,'LineWidth',2,Color=ipe.darkred)
hold off;
legend('20C','10C','0C','-10C','-20C','fit','Location','northwest')
xlabel('f_{master}')
ylabel('Dynamic modulus master')

out = false;
end

The error message then looks like this:

Answer 1

[Ishu]

Hi Ondrej,

I understand that you are getting error while using "solve()" function. This is because your function "plotprogress(...)" is not valid as in this function you are plotting at 3 different axes in one iteration and this is not supported by "plotFcn" . And as you said that when you replace "subplot(...)" with "gca" it does not give any error, this is because "gca" uses only one set of axes to plot.

So what use can do is:

1. You can replace "PlotFcn" with "OutputFcn"

lsqoptions.OutputFcn = @(x, optimValues, flag) plotprogress(x, optimValues, flag, ipe, plotting, ord, freq, stor_mod, loss_mod, dyn_mod, E);

2. You can create different functions for different plots, and provide all these functions in "lsqoptions.PlotFcn". This will take more time to complete.

lsqoptions.PlotFcn = {@(x, optimValues, flag) plotprogress_storage(..), @(x, optimValues, flag) plotprogress_loss(..), @(x, optimValues, flag) plotprogress_dynamic(..)};

function out  = plotprogress_storage(..)
subplot(3,1,1)

for 
    % your code
end

% your code 
end

function out  = plotprogress_loss(..)
subplot(3,1,2)

for 
    % your code
end

% your code 

end

function out  = plotprogress_dynamic(..)
subplot(3,1,3)

for 
    % your code
end

% your code 

end

For more information on "lsqnonlin" you can refer to this documentation:

https://www.mathworks.com/help/optim/ug/lsqnonlin.html

In the above documentation you can also check for different "options" available in "lsqnonlin".

For Output function and plot function you can refer this:

https://www.mathworks.com/help/optim/ug/output-function.html

Hope it helps!