Genetic Algorithm output is different than manual calculation

below is a code to minimize a function y using genetic algorithm.

function P= decay_one(k,n)
y=1000000*exp(0.000005*x).*sin(2*pi*5000*x.^2);
end
lb=0;
ub=100;
n=1;
[x,fval]=ga(@decay_one,n,[],[],[],[],lb,ub);

-----------------------------------------------

GA output

x = 78.4000, fval = -9.9894e+05

-----------------------------------

But If i plot the same in matlab manually using below code, I get response as attached below.

t=1:1:100;
S=1000000*exp(0.000005*t).*sin(2*pi*5000*t.^2);
plot(t,S)

------------------------------------------

In the plot the minimum value of Y seems to be somehwere btn -0.02 to -0.03, but GA minimum comes out to be fval = -9.9894e+05.

PLEASE HELP ME UNDERSTAND THE VARIATION IN MANUAL PLOT AND GA OUTPUT.

Answer 1

[Sam Chak]

Hi @Kalyan Dash

Look at the angular frequency of the sinusoidal component, you will see it is a high-frequency signal.

In your code, the sampling rate is very coarse t = 1:1:100;. Thus, it cannot capture all details in that spectrum.

If the step size is much smaller, then you can see the plot.

x = 99.9:0.000001:100;
y = 1000000*exp(0.000005*x).*sin(2*pi*5000*x.^2);
plot(x, y)

Answer 2

[Alan Weiss]

The issue is that your objective function varies extremely quickly, with about 10000 squiggles per unit near x = 1. The exponential term does verly little except to ensure that any optimum you find is more likely to be near the upper end of your range than the lower end. Near your upper bound of 100 the function is oscillating so rapidly that you cannot plot it accurately on the scale you have chosen. Try this:

t = linspace(99,99.00001); % Equally spaced points from 99 to 99.00001
plot(t,1000000*exp(0.000005.*t).*sin(2*pi*5000*t.^2));

Do you see how your exponential term does nothing, and the function oscillates so rapidly that you are essentially sampling random noise?

Alan Weiss

MATLAB mathematical toolbox documentation

Answer 3

[Walter Roberson]

S=1000000*exp(0.000005*t).*sin(2*pi*5000*t.^2);

Consider sin(2*pi*5000*t): that would be a sine wave with frequency 5000 cycles per 2*pi, so with sin(2*pi*5000*t) you would expect 5000 cycles approximately every 6.3, or approximately 796 cycles per unit. But your t is 1:100 so you are sampling less than 1/796 th of the peaks per unit.

But you do not have sin(2*pi*5000*t) you have sin(2*pi*5000*t^2) . Does that increase or decrease the frequency? Clearly it increases the frequency -- by 75-ish your sine frequency is over 2 megahertz.

Where is the true minima?

Well, exp(0.000005*t) is exponential increasing in t, and sin() of real values stays in the range -1 to +1 so for any t that gives a positive S, a value of t that leads to exactly one cycle later in the sine wave would give the same sin() value but would have increasing exp(0.000005*t) (because larger t) and so would lead to a larger value of S. Likewise, for any given t that gives a positive S, there is a slightly larger value of t1 that is going to lead to the sin() being -1, and then there would be a t2 slightly larger than that which would lead to the same point on the sin() cycle but with t2>t1 then exp(0.000005*t2) > exp(0.000005*t1) and it follows that S for t2 would be more negative than the S for t1.

We thus see that no matter what minima we find for S, there is a time nearby that will give more of a minima. Therefore the true minima is arbitrarily small.

(In this case you would not say that the minima is -infinity at t = infinity because sin(infinity) is not defined.)

Answer 4

[Kalyan Dash]

Thanks for your swift response @Sam Chak, @Alan Weiss and @Walter Roberson. All of you have rightly pointed out the mistake I was doing. The range of 0:100 is too large for the high frequency signal. Therefore, the manual plot was unable to capture all mimima pooints accurately. However, GA was giving close to correct answer.

[Sam Chak]

Hi @Kalyan Dash

Theoretically, the minimum point should locate at the last trough of the search interval.

% Find minimum point
fun = @decay_one;
lb = 0;     % lower bound
ub = 3;     % upper bound
n  = 1;
[xsol, fval] = ga(fun, n, [], [], [], [], lb, ub)

Optimization terminated: average change in the fitness value less than options.FunctionTolerance.

xsol = 2.9588

fval = -19.2671

% Plot function
x = linspace(0, 3, 3001);
y = 1*exp(1*x).*sin(2*pi*1*x.^2);
plot(x, y), grid on, hold on
plot(xsol, fval, 'o', 'linewidth', 2, 'MarkerSize', 10)
xlabel('x'), ylabel('y')

function P = decay_one(x)
    P = 1*exp(1*x).*sin(2*pi*1*x.^2);
end