## Problem with including ode45 in paretosearch

### henry kagey (view profile)

on 30 Jan 2019
Latest activity Edited by henry kagey

### henry kagey (view profile)

on 31 Jan 2019
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 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
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

### Alan Weiss (view profile)

on 31 Jan 2019
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

### henry kagey (view profile)

on 31 Jan 2019
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.
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.