Problem with including ode45 in paretosearch

Hello all,

I had a perectly functioning paretoserach script that I tried to improve by simultaneously solving a differential equations using ode45 inside the solver. I tested the ode45 code separately, and it does exactly what I expected it to do. I also tested it inside the longer algorithm over a single run (before applying the paretosearch) anbd again it worked flawlessly. After adding this into the paretosearch script the code runs, but the number of pareto solutions is small, and hardly spaced apart at all. I also used (an already functioning) gamultiobj script with the same ode45 inclusion, and while I get a number of solutions using this code, they are only infinitessimally spaced apart.

i should also mention that the previous run time of the paretosearch was on the order of 20 seconds. After including ode45 the time to solve is on the order of 20 minutes. The original function count is on the order of 2000, and after ode45 is added the function count is 15000 (and appears to be maxxing out).

Below is an image of the paretosearch code functioning before adding in ode45, and then the result after ode45 is included:

After searching about some, I put the ode45 call inside of a (nested) function within the paretosearch function, but this did not change anything in the solution. I do not undertand what could be causing the problem. I have included a pseudocode below in order to clarify what I am asking about specifically.

%% Paretosearch and ode45 pseudocode


%% Section 0

% Initial definitions

% start timer
tic

%% Section 0.1: definitions 

%define vector c1(n)

%define vector c2(n)

%% Section 1: definition of parameters for paretosearch

% parameters of interest
% x(1) = tilt (deg)
% x(2) = diameter (m)
% x(3) = depth (m)

lb = [1, 0.005, 0.02];
ub = [89, 0.015, 0.1];
% lower and upper bounds on design variables,

fun = @(x)objval(x, c1, c2);  %objective function handle

%% Section 1.1: paretosearch
npts = 200;
opts_ps.ParetoSetSize = npts;

%options for paretosearch
opts_ps = optimoptions('paretosearch','Display','off','PlotFcn','psplotparetof');

% definition of contstraints and number of variables

A = []; b = [];
Aeq = []; beq = [];
numberOfVariables = 3;

% function paretosearch call

% Note: calls 'fun' as one of the variables it passes
[x_ps,fval_ps,exitflag,psoutput1] = paretosearch(fun,numberOfVariables,A,b,Aeq,beq,lb,ub,[],opts_ps);

disp("NSGA-II Total Paretosearch Function Count: " + psoutput1.funccount);


% find optimal points to start pareto solution
x0 = zeros(2,3);
x0f = (lb + ub)/2;
opts_fmc = optimoptions('fmincon','Display','off','MaxFunctionEvaluations',1e4);
x0(1,:) = fmincon(@(x)pickindex(x,1,c1,c2),x0f,[],[],[],[],lb,ub,[],opts_fmc);
x0(2,:) = fmincon(@(x)pickindex(x,2,c1,c2),x0f,[],[],[],[],lb,ub,[],opts_fmc);

opts_ps.InitialPoints = x0;
opts_ps.PlotFcn = [];
[x_psx0,fval_ps1x0,exitflag,psoutput1x0] = paretosearch(fun,numberOfVariables,A,b,Aeq,beq,lb,ub,[],opts_ps);
disp("Total Initial Opt pt. Paretosearch Function Count: " + psoutput1x0.funccount);

%resolve with optimal start points and weighted solver
Fmax = max(fval_ps1x0);
nobj = numel(Fmax);
Fmin = min(fval_ps1x0);
w = sum((Fmax - fval_ps1x0)./(1 + Fmax - Fmin),2);
p = w.*((Fmax - fval_ps1x0)./(1 + Fmax - Fmin));
xnew = zeros(size(x_psx0));
nsol = size(xnew,1);
fvalnew = zeros(nsol,nobj);
opts_fg = optimoptions('fgoalattain','Display','off');
nfv = 0;
for ii = 1:nsol
    [xnew(ii,:),fvalnew(ii,:),~,~,output] = fgoalattain(fun,x_psx0(ii,:),fval_ps1x0(ii,:),p(ii,:),...
        A,b,[],[],lb,ub,[],opts_fg);
    nfv = nfv + output.funcCount;
end
disp("fgoalattain Function Count: " + nfv)


% Plotting/comparing solutions 
figure
fps = sortrows(fval_ps,1,'ascend');
plot(fps(:,1),fps(:,2),'g-')
hold on
fps2 = sortrows(fval_ps1x0,1,'ascend');
plot(fps2(:,1),fps2(:,2),'r-')

[fps3,indexp] = sortrows(fvalnew,1,'ascend');
plot(fps3(:,1),fps3(:,2),'k.-')

legend('fpval_ps','fval_ps1x0','fvalnew')
%'paretoserach','paretoserach w/inital pts', 'paretosearch w/hybrid')

xlabel('objective(f1)')
ylabel('objective(f2)')
title ('Paretosearch algorithm for solution output')

toc;
time = toc

%% function pick index for inital optimal points
function z = pickindex(x,k,c1,c2)
    z = objval(x, c1, c2); % evaluate both objectives
    z = z(k); % return objective k
end

    
%% Section 2: objective function defintion

function F = objval(x, c1, c2)
% here I pass c1 and c2 along with x into the evaluation fo the objective function
%nargin(objval)
%% Section 2.1: physical parameters

% here I define a number of parameters which depend on the parameters of interest

%% Section 2.2 Forcing function data 
% forcing function acts on physical object

elez = zeros;
aziz = zeros;
[ele,azi] = data1(elez,aziz);
% ele = elevation azi = azimuth

% length of data vector
n = length(ele);
delta = const;

[alpha,gamma] = datatovector1(ele,azi,n,delta);
%alpha and gamma are used to define theta in the next function

sigma = degtorad(x(1));
% sigma is the angle to the horizontal plane

[theta] = incidentangle1(alpha,gamma,n,sigma);


%% Section 2.3: Source values at specific times

%define more constants and vectors


%% Section 2.4: fraction of forcing function due to geometry
% Angle relations for removal

% this section relates the angle of the incident radiation on the panel
% to the angle in the cylinder to determine the fraction of cylinder
% surface area illuminated
for i = 1:l
    thetagl(i) = asin( 1/1.5 * sin(theta(i)) );
    if thetagl(i) < atan(dcyl/(dcyl+lcyl));
        Acyl(i) = pi*dcyl*lcyl/2;
    else
        Acyl(i) = 4 * (dcyl/2)^2 * cot(thetagl(i));
    end
    Afrac(i) = Acyl(i)/(pi*dcyl*lcyl); 
    Atot(i) = Afrac(i)*Areacyl*numcyl;
end

I_o = 1000;

for i = 1:l;
    ang = (round(radtodeg(theta(i))));%getting a count number
    ta = c2(ang);
    Ifracb(i) = I_o*ta/187.4 * ( sin(theta(i)) / cot(thetagl(i)) ) * (pi/4);
    % 
    Ifracph(i) = ( (I_o*cos(theta(i))*ta)+100 )/446;
    % 
end



%% Section 2.5: Losses 

% define more constants

%% Section 2.6: physical study


delTo = Quo*600/(const); %seconds to 10 minute time intervals
T(1) = c1(1)+delTo;

Qm = const*x(3);
Vpanel = some value; %f(x(i))

Ctot(1) = 1;

kphconst = const;

count = 0;

%initialize output objective functions
Vremtot = 0;
Egaintot = 0;

for i = 1:l-1;
    %initialize phsyical values
    
    %...multiple calulations...
    
    T_o = T(i)+delT_o;
    Tavg(i) = (T_o+T(i))/2;
    
    %iteration for more accurate values
    
    %...multiple calulations...
    
    T(i+1) = T(i)+delT(i);
    
    
    % begin primary study
    
    % define constants
    
    kvb = some value; %f(x(i))
    kvd = some value; %f(x(i))
    kph = some value; %f(x(i))
   
    ktot = (kvb + kvd + kph);
    
    
   %% ode45 to solve simultaneous ODEs
   
    C_old = Ctot(i);
    
    f = @(t,y)[ (Qm/Vpanel) * ( y(1)-y(2) ) - ktot*y(2) ; (Qm/Vres) * ( y(1)-y(2) ) ];
    tspan = [0 10];
    y1_o = C_old;
    y2_o = C_old;
    y0 = [y1_o , y2_o];
    [Time Y] = ode45(f,tspan,y0);
    
    last = length(Y(:,2));
    Ctot(i+1) = Y(last,2);
    
    clear Y;
    
    %This section for summing  objective functions
    if Ctot(i+1)<=10^-3;
        Ctot(i+1) = 1;
        Vremtot = Vremtot + Vtot;
        Egaintot = (T(i+1)-c1(i))*(cptot*Mtot) + Egaintot;
        T(i+1) = c1(i)(i+1);
        count = count+1;
    else
        continue
    end
    
end

% objective functions
Vremtot = Vremtot + (1-Ctot(l))*Vtot;
Egaintot = Egaintot + (T(l)-c1(l))*(cptot*Mtot);

F = [-Vremtot, -Egaintot];
% this is the return of the objective function evaluation 'objval'


%% Functions within objval function

% functions for data, included inside objective function evaluation
% because the tilt of the panel is one of the variables of interest (x(1)).

%% function data1
    function [ele,azi] = data1(elez,aziz)
        
        ele = [some vector];
        
        azi = [some vecotr];
        
    end


%% function datatovector1

    function [alpha,gamma] = datatovector1(ele,azi,n,delta);
        
        for i = 1:n;
            
            elea(i) = 90 - ele(i);
            
            if elea(i) <= 0;
                continue;
            end
            
            if abs(delta - azi(i)) >= 90;
                continue
            end
            
            if elea(i) >= 90;
                elea(i) = 90;
            end
            
            
            alpha(i) = degtorad(elea(i));
            gamma(i) = degtorad(azi(i) - delta);
        end
        
    end

%% function incidentangle1

    function [theta] = incidentangle1(alpha,gamma,n,sigma)
        
        for i = 1:n-1;
            
            costheta(i) = sin(alpha(i))*cos(sigma) + cos(gamma(i))*cos(alpha(i))*sin(sigma);
            
            theta(i) = acos(costheta(i));
            
        end
    end




end
% end of functions within objval function

Any help or assistance with this problem would be much appreciated. I have been beating my head against this for a week now with no results.

Thank you

[Alan Weiss]

Can you please explain what your ode45 call is supposed to do, how does it relate to your original multiobjective optimization? I'd like some context for it so I can understand why you think that an ode45 call should improve the behavior of the multiobjective solvers. It is clear that the call is having a negative effect, but I don't understand how it relates at all.

Alan Weiss

MATLAB mathematical toolbox documentation

[henry kagey]

Hi Alan,

The ode45 part of the script calculates concentrations at two points in a reactor. The paretosearch part of the script it is embedded in is investigating geometric parameters related to reactor performance.

I was thinking about this and wonder if the function handle in the ode45 call:

f = @(t,y)[ (Qm/Vpanel) * ( y(1)-y(2) ) - ktot*y(2) ; (Qm/Vres) * ( y(1)-y(2) ) ];

is disrupting the anonymous function handle call in used for the pareto search:

fun = @(x)objval(x, c1, c2);

I'm not sure why this might have this affect, but perhaps it is? It appears the size of the solution vector, or the pareto set size, is being constrained after including the ode45 call.