Asked by henry kagey
on 30 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 = 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

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Alan Weiss (view profile)

Direct link to this comment:https://de.mathworks.com/matlabcentral/answers/442361-problem-with-including-ode45-in-paretosearch#comment_666384

## henry kagey (view profile)

Direct link to this comment:https://de.mathworks.com/matlabcentral/answers/442361-problem-with-including-ode45-in-paretosearch#comment_666394

Sign in to comment.