How to do optimization for result calculated by using experimental data?

Hi All,

I have run a measurement for 28 fabricated samples based on few geometrical parameter (x(1), x(2), x(3) and x(4)). Each sample has a different combination of geometrical parameters. From the measurement result, i got list of 'alpha' and 'xi' (that varied with frequency) for each sample. Based on the 'alpha' and 'xi' values, I calculate the output, y (versus frequency) for each sample. With the data I have, I would like to find the optimized parameter for x(1), x(2), x(3) and x(4) when the output is maximized to 1 (y=1) at frequency 1600Hz.

I tried a few method but unfortunately it fails to achieve my goal.

I really appreaciate if someone can help me to solve this problem.

I have attached the files that I have for your kind reference. Kindly, please advice me the next step to achieve my goal.

Thank you.

Regards,

Nur Arafah

alpha = load ('alpha.txt'); %% alpha
xi = load ('xi.txt');  %% xi
c = 345.2;
n = 1.83e-5;
den = 1.189;    
f = 500:1.5625:4000; %% frequency

%%  for samples with different x(1)
x(1)=700e-6; 
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,1);
    E(i) = xi(i,1);
    x(2) = 1000e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2;  %% output for sample #1
end
y1=y;

x(1)=750e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,2);
    E(i) = xi(i,2);
    x(2) = 1000e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #2
end
y2=y;

x(1)=800e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,3);
    E(i) = xi(i,3);
    x(2) = 1000e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #3
end
y3=y;

x(1)=850e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,4);
    E(i) = xi(i,4);
    x(2) = 1000e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #4
    end
y4=y;

x(1)=900e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,5);
    E(i) = xi(i,5);
    x(2) = 1000e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #5
end
y5=y;

x(1)=950e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,6);
    E(i) = xi(i,6);
    x(2) = 1000e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #6
end
y6=y;

x(1)=1000e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,7);
    E(i) = xi(i,7);
    x(2) = 1000e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #7
end
y7=y;

%% for samples with different x(2)
x(2)=700e-6; 
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,8);
    E(i) = xi(i,8);
    x(1) = 700e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #8
end
y8=y;

x(2)=800e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,9);
    E(i) = xi(i,9);
    x(1) = 700e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #9
end
y9=y;

x(2)=900e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,10);
    E(i) = xi(i,10);
    x(1) = 700e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #10
end
y10=y;

x(2)=1000e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,11);
    E(i) = xi(i,11);
    x(1) = 700e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #10
end
y11=y;

x(2)=1100e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,12);
    E(i) = xi(i,12);
    x(1) = 700e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #12
end
y12=y;

x(2)=1200e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,13);
    E(i) = xi(i,13);
    x(1) = 700e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #13
end
y13=y;

x(2)=1300e-6;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,14);
    E(i) = xi(i,14);
    x(1) = 700e-6
    x(3) = 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #14
    end
y14=y;

%%  for samples with different x(3)
x(3)=0.004; 
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,15);
    E(i) = xi(i,15);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #15
    end
y15=y;

x(3)=0.005;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,16);
    E(i) = xi(i,16);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #16
    end
y16=y;

x(3)=0.006;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,17);
    E(i) = xi(i,17);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #17
   end
y17=y;

x(3)=0.007;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,18);
    E(i) = xi(i,18);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #18
    end
y18=y;

x(3)=0.008;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,19);
    E(i) = xi(i,19);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #19
    end
y19=y;

x(3)=0.009;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,20);
    E(i) = xi(i,20);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #20
    end
y20=y;

x(3)=0.01;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,21);
    E(i) = xi(i,21);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #21
end
y21=y;

%%  for samples with different x(4)
x(4)=0.001; 
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,22);
    E(i) = xi(i,22);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #22
end
y22=y;

x(4)=0.002;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,23);
    E(i) = xi(i,23);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #23
end
y23=y;

x(4)=0.003;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,24);
    E(i) = xi(i,24);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #24
end
y24=y;

x(4)=0.004;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,25);
    E(i) = xi(i,25);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #25
end
y25=y;

x(4)=0.005;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,26);
    E(i) = xi(i,26);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #26
end
y26=y;

x(4)=0.006;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,27);
    E(i) = xi(i,27);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #27
end
y27=y;

x(4)=0.007;
for i = 1:length(f)
    w(i) = 2.*pi.*f(i);
    a(i) = alpha(i,28);
    E(i) = xi(i,28);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    y(i) = 1 - abs(RF(i)).^2; %% output for sample #28
end
y28=y;

%% all the output result from sample #1 to #28
f=[f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f;f]
y=[y1;y2;y3;y4;y5;y6;y7;y8;y9;y10;y11;y12;y13;y14;y15;y16;y17;y18;y19;y20;y21;y22;y23;y24;y25;y26;y27;y28]
plot(f,y);

[Mario Malic]

If I understood your code correctly, you have vector for each sample that contains output values y for the frequencies in range [500, 4000].

You just need to find the value of i_freq that corresponds to frequency of 1600 Hz and find the maximum value of y(i_freq) for all the vector samples.

Since the code is quite repetitive, the only varibles that change are x(1) to x(4), write a function that you call for every sample

% You send a vector f and x variables, and receive y as a output vector and x as a vector of values
[y, x] = outputcalc(f, x(1), x(2), x(3), x(4)) 

Then you can save it in a array that each sample represents a row of a matrix. Then your problem would be just finding max value of a specific column in y, and with the index of it, you call the particular row of x array.

[Jafar Nur Arafah]

Hi Mario Malic,

Thanks for the suggestion.

Based on what I understood in your comment, you suggest to write a function 'outputcalc' for every sample based on the save file of y, f, x(1), x(2), x(3) and x(4). Then call the function to represent each sample based on representative row of matrix. However, when I tried to do so, the error came out and mentioned "Unrecognized function or variable 'outputcalc'."

FYI, at the mean time. I have try to tune the code that I have and find the optimization by using GA Optimization tool box. I focused only on the frequency 1600Hz. I maximize the output/objective by adding '-' in the code.

And finally I success to get the optimized parameter that I want. But, I found out that the predicted output given by Matlab is too far from the target (y=1), which is -0.4410947. Where if I reverse the result it will be 0.4410947.

What can add in the code, so that the objective function is maximized to certain value? Example like my target, I want the output to maximize the predicted output value to 1.

function [y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11,y12,y13,y14,y15,y16,y17,y18,y19,y20,y21,y22,y23,y24,y25,y26,y27,y28]=file2_3(x)


alpha = load ('alpha.txt');
xi = load ('xi.txt');
c = 345.2;
n = 1.83e-5;
den = 1.189;    
f1=1600;

%% d
x(1)=700e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,1);
    E(i) = xi(705,1);
    x(2)= 1000e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #1
    y(i)= -B(i); %% maximize the output/objective
    
end
y1=y;

x(1)=750e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,2);
    E(i) = xi(705,2);
    x(2)= 1000e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #2
    y(i)= -B(i); %% maximize the output/objective
     
end
y2=y;

x(1)=800e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,3);
    E(i) = xi(705,3);
    x(2)= 1000e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #3
    y(i)= -B(i); %% maximize the output/objective
    
end
y3=y;

x(1)=850e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,4);
    E(i) = xi(705,4);
    x(2)= 1000e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #4
    y(i)= -B(i); %% maximize the output/objective
    
    end
y4=y;

x(1)=900e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,5);
    E(i) = xi(705,5);
    x(2)= 1000e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #5
    y(i)= -B(i); %% maximize the output/objective
    
end
y5=y;

x(1)=950e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,6);
    E(i) = xi(705,6);
    x(2)= 1000e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #6
    y(i)= -B(i); %% maximize the output/objective
    
end
y6=y;

x(1)=1000e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,7);
    E(i) = xi(705,7);
    x(2)= 1000e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #7
    y(i)= -B(i); %% maximize the output/objective
    
end
y7=y;

%% t
x(2)=700e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,8);
    E(i) = xi(705,8);
    x(1)= 700e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #8
    y(i)= -B(i); %% maximize the output/objective
    
end
y8=y;

x(2)=800e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,9);
    E(i) = xi(705,9);
    x(1)= 700e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #9
    y(i)= -B(i); %% maximize the output/objective
    
end
y9=y;

x(2)=900e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,10);
    E(i) = xi(705,10);
    x(1)= 700e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #10
    y(i)= -B(i); %% maximize the output/objective
    
end
y10=y;

x(2)=1000e-6
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,11);
    E(i) = xi(705,11);
    x(1)= 700e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #11
    y(i)= -B(i); %% maximize the output/objective
    
end
y11=y;

x(2)=1100e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,12);
    E(i) = xi(705,12);
    x(1)= 700e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #12
    y(i)= -B(i); %% maximize the output/objective
    
end
y12=y;

x(2)=1200e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,13);
    E(i) = xi(705,13);
    x(1)= 700e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #13
    y(i)= -B(i); %% maximize the output/objective
    
end

y13=y;

x(2)=1300e-6;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,14);
    E(i) = xi(705,14);
    x(1)= 700e-6
    x(3)= 0.00615
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #14
    y(i)= -B(i); %% maximize the output/objective
    
    end
y14=y;

%% p
x(3)=0.004;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,15);
    E(i) = xi(705,15);
    x(1)= 700e-6
    x(2)= 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #15
    y(i)= -B(i); %% maximize the output/objective
    
    end
y15=y;

x(3)=0.005
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,16);
    E(i) = xi(705,16);
    x(1)= 700e-6
    x(2)= 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #16
    y(i)= -B(i); %% maximize the output/objective
    
    end
y16=y;

x(3)=0.006;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,17);
    E(i) = xi(705,17);
    x(1)= 700e-6
    x(2)= 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #17
    y(i)= -B(i); %% maximize the output/objective
    
   end
y17=y;

x(3)=0.007;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,18);
    E(i) = xi(705,18);
    x(1)= 700e-6
    x(2)= 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #18
    y(i)= -B(i); %% maximize the output/objective
    
    end
y18=y;

x(3)=0.008;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,19);
    E(i) = xi(705,19);
    x(1)= 700e-6
    x(2)= 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #19
    y(i)= -B(i); %% maximize the output/objective
    
    end
y19=y;

x(3)=0.009;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,20);
    E(i) = xi(705,20);
    x(1)= 700e-6
    x(2)= 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #20
    y(i)= -B(i); %% maximize the output/objective
    
    end
y20=y;

x(3)=0.01;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,21);
    E(i) = xi(705,21);
    x(1)= 700e-6
    x(2)= 1000e-6
    x(4) = 0.0042
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #21
    y(i)= -B(i); %% maximize the output/objective
    
end
y21=y;

%% D
x(4)=0.001; 
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,22);
    E(i) = xi(705,22);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #22
    y(i)= -B(i); %% maximize the output/objective
    
end
y22=y;

x(4)=0.002;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,23);
    E(i) = xi(705,23);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #23
    y(i)= -B(i); %% maximize the output/objective
    
end
y23=y;

x(4)=0.003;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,24);
    E(i) = xi(705,24);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #24
    y(i)= -B(i); %% maximize the output/objective
    
end
y24=y;

x(4)=0.004;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,25);
    E(i) = xi(705,25);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #25
    y(i)= -B(i); %% maximize the output/objective
    
end
y25=y;

x(4)=0.005;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,26);
    E(i) = xi(705,26);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #26
    y(i)= -B(i); %% maximize the output/objective
    
end
y26=y;

x(4)=0.006;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,27);
    E(i) = xi(705,27);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #27
    y(i)= -B(i); %% maximize the output/objective
    
end
y27=y;

x(4)=0.007;
for i = 1:length(f1)
    w(i) = 2.*pi.*f1;
    a(i) = alpha(705,28);
    E(i) = xi(705,28);
    x(1) = 700e-6
    x(2) = 1000e-6
    x(3) = 0.00615
    k(i) = x(1).*sqrt((w(i).*den)./(4.*n));
    Rs(i) = (1./2).*sqrt(2.*n.*den.*w(i));
    r1(i) = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k(i).^2)./32)))+((2.*a(i).*Rs(i))./(x(3).*den.*c));
    m1(i) = ((((w(i).*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k(i).^2)./2)))))+(((w(i).^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w(i))./(3.*pi.*x(3).*c)).*E(i))-cot((w(i).*x(4))./c);
    z(i) = r1(i)+(j.*m1(i));
    Ref(i) = (z(i)-1)./(z(i)+1);
    Ref_re(i) = real(Ref(i));
    Ref_im(i) = imag(Ref(i));
    RF(i) = sqrt((Ref_re(i).^2)+(Ref_im(i).^2));
    B(i) = 1 - abs(RF(i)).^2; %%calculated output for sample #28
    y(i)= -B(i); %% maximize the output/objective
    
end
y28=y;


y=[y1;y2;y3;y4;y5;y6;y7;y8;y9;y10;y11;y12;y13;y14;y15;y16;y17;y18;y19;y20;y21;y22;y23;y24;y25;y26;y27;y28];
end

[Mario Malic]

• Please see the attached files. There's no need for this simple script to be this long.
• I hope the code represents what you want to do, in the other words, check if you get the same outputs
• For your and our sanity, work with arrays.

x = [700e-6, 1000e-6, 0.00615, 0.0042]; 

y = Sample_Calculation (x, f, alpha, xi, c, n, den, Sample_num);

EDIT 1: (Question is: what is your output? y is a (ixi) matrix and you mentioned how you want to maximize y on 1600 Hz. Most definitely is that you have expected a y as a vector so it means that you didn't write your code well.) I didn't clear variables, error on my side. Will update

EDIT 2: The second question is, do you need to use the optimization algorithm? Are you looking for minimum of y for given x(1), x(2), x(3), x(4) variables? Didn't read well either.

[Mario Malic]

Here should be something good for you, at least to move further with your code.

I assume x variables need to be positively defined, so you need to include that in your constraints.

Can't upload more files, as I don't know how to delete the ones which are uploaded before, that's why it's written like this.

clc;
clear vars;
alpha = load ('alpha.txt');
xi = load ('xi.txt');
c = 345.2;
n = 1.83e-5;
den = 1.189;    
f = 500:1.5625:4000;
Sample_Num = 1; % This doesn't change, unnecessary probably
Desired_Freq = find (f == 1600);

%% Calculating objective value for each sample
options = optimoptions('ga','TolFun', 1e-3, 'PlotFcn', @gaplotbestf);
x = [x(1), x(2), x(3), x(4)];
[X_Opt, fval, exitflag,output] = ga(@(x)Sample_Calculation(x, f, alpha, xi, c, n, den, Sample_Num, Desired_Freq), 4, options)
% y = Sample_Calculation (X_Opt, f, alpha, xi, c, n, den, Sample_Num, Desired_Freq)

Sample_Calculation function

function y = Sample_Calculation (x, f, alpha, xi, c, n, den, Sample_Num, Desired_Freq)
w = 2.*pi.*f;
a = alpha(:,Sample_Num)';
E = xi(:,Sample_Num)';
k = x(1).*sqrt((w.*den)./(4.*n));
Rs = (1./2).*sqrt(2.*n.*den.*w);
r1 = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k.^2)./32)))+((2.*a.*Rs)./(x(3).*den.*c));
m1 = ((((w.*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k.^2)./2)))))+(((w.^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w)./(3.*pi.*x(3).*c)).*E)-cot((w.*x(4))./c);
z = r1+(j.*m1);
Ref = (z-1)./(z+1);
Ref_re = real(Ref);
Ref_im = imag(Ref);
RF = sqrt((Ref_re.^2)+(Ref_im.^2));
yn= 1 - abs(RF).^2;  %% output for sample #1
y = 1 + yn(1,Desired_Freq); % Objective function, please CHECK this
% plot (f,yn); % only displays last iteration 
end

[Jafar Nur Arafah]

Hi Mario Malic,

Thanks for the support.

Below are my reply to your answer above:

Q: I hope the code represents what you want to do, in the other words, check if you get the same outputs

A: Firstly, I have cross check the output (y vs frequency) from the Jafar.m file. The output does not match with the output from the first code I upload (you can refer to file1.m attached). The result from Jafar.m is a matrix of 2241x2241, while the result from file1.m is a matrix of 28x2241. Where 28 represent the 28 samples and 2241 represent the result (y) from 500-4000Hz.

EDIT 1: (Question is: what is your output? y is a (ixi) matrix and you mentioned how you want to maximize y on 1600 Hz. Most definitely is that you have expected a y as a vector so it means that you didn't write your code well.) I didn't clear variables, error on my side. Will update

A: My output is y (varied across frequency), from file1.m you can see that the frequency range is in between 500-4000Hz. However, as my goal is to maximize the output (y) at 1600Hz. You can see that re-write code in my second post (you can refer to file2_3.m attached) and focused only for output (y) at 1600Hz for each sample and do optimization from the result.

EDIT 2: The second question is, do you need to use the optimization algorithm? Are you looking for minimum of y for given x(1), x(2), x(3), x(4) variables? Didn't read well either.

A: Based on the output of 28 sample that I have, I would like to find the optimized x(1), x(2), x(3) and (4) value that will give maximum of output (y) equal to 1 (at frequency 1600Hz).

I assume x variables need to be positively defined, so you need to include that in your constraints.

As you know each sample have different x(1), x(2), x(3) and x(4) value (you can refer to Table.JPG). From your comment you mentioned that I need to defined the x variable in the code by using constraint method. However, I think it is not proper to define the x variable in constraint. This is because, the x variable is connected to the alpha and xi values which is specified for each sample. Define the x variables in constraint method might lead to different output for each sample and might effected the optimization result.

Thanks alot for your advices and suggestions. It really help the Matlab's amateur like me. =)

Regards,

Nur Arafah

[Mario Malic]

Here's code which is not working because y is not equal to the same variable from your code. I can't spend time troubleshooting where did things go wrong so, you have to figure it out yourself.

Also, reading your post history, I doubt this is what you actually wanted. Since this only gives you a maximum of your experiments, not the actual maximum.

format long
clc;
clearvars;

alpha = (load ('alpha.txt'))';
xi = (load ('xi.txt'))';
c = 345.2;
n = 1.83e-5;
den = 1.189;    
f = 500:1.5625:4000;
Sample_Num = 28;
Desired_Freq = find (f == 1600);

% Defining X matrix
x = zeros(28,4);
x(1:7, 1) = [700e-6:50e-6:1000e-6];
x(8:28,1) = 700e-6;
x(1:7,2) = 1000e-6;
x(8:14, 2) = [700e-6:100e-6:1300e-6];
x(15:28,2) = 1000e-6;
x(1:14, 3) = 0.00615;
x(15:21, 3) = [0.004:0.001:0.010];
x(22:28,3) = 0.00615;
x(1:21, 4) = 0.0042;
x(22:28, 4) = [0.001:0.001:0.007];

y = zeros (Sample_Num, length(f));
for ii = 1 : 1 : Sample_Num
y(ii,:) = Sample_Calculation (x(ii,:), ii, f, alpha, xi, c, n, den);
end
[Max_Y, Index] = max(y(:,Desired_Freq))
X_Opt = x(Index,:)

Function

function y = Sample_Calculation (x, ii, f, alpha, xi, c, n, den)
w = 2.*pi.*f;
a = alpha(ii,:);
E = xi(ii,:);
k = x(1).*sqrt((w.*den)./(4.*n));
Rs = (1./2).*sqrt(2.*n.*den.*w);
r1 = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k.^2)./32)))+((2.*a.*Rs)./(x(3).*den.*c));
m1 = ((((w.*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k.^2)./2)))))+(((w.^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w)./(3.*pi.*x(3).*c)).*E)-cot((w.*x(4))./c);
z = r1+(j.*m1);
Ref = (z-1)./(z+1);
Ref_re = real(Ref);
Ref_im = imag(Ref);
RF = sqrt((Ref_re.^2)+(Ref_im.^2));
y = 1 - abs(RF).^2;  %% output for sample #1
end

[Jafar Nur Arafah]

Hi Mario Malic,

Thanks for your time. Really appreciate it.

The code works well. But as you said it find the maximum of y based on the data plug in.

However, my goal to use all the y data I have and do the optimization (y=1) at frequency. When I said use all the data, meaning that I want to combine all the data but it is not sum(y).

So, what I do is, first I focus the output only at 1600Hz (instead of 500-4000Hz). Where in the alpha and xi file, I remove other frequency value and maintain the result for 1600Hz only.

%% Sample_Calculation_3
function y = Sample_Calculation_3 (x, ii, f, alpha, xi, c, n, den)
w = 2.*pi.*f;
a = alpha(ii,:);
E = xi(ii,:);
k = x(1).*sqrt((w.*den)./(4.*n));
Rs = (1./2).*sqrt(2.*n.*den.*w);
r1 = (((32.*n.*x(2))./(x(3).*den.*c.*(x(1).^2))).*sqrt(1+((k.^2)./32)))+((2.*a.*Rs)./(x(3).*den.*c));
m1 = ((((w.*x(2))./(x(3).*c)).*(1+(1./sqrt(9+((k.^2)./2)))))+(((w.^2).*(x(1).^2))./(8.*x(3).*(c.^2)))+((8.*x(1).*w)./(3.*pi.*x(3).*c)).*E)-cot((w.*x(4))./c);
z = r1+(j.*m1);
Ref = (z-1)./(z+1);
Ref_re = real(Ref);
Ref_im = imag(Ref);
RF = sqrt((Ref_re.^2)+(Ref_im.^2));
y = 1 - abs(RF).^2;  %% output for sample #1
end

Then, I edit the code u give me and make it become the objective function.

%% Objective function
function y = Fifth_edit(x)
format long
clc;
clearvars;

alpha = (load ('alpha.txt'))';
xi = (load ('xi.txt'))';
c = 345.2;
n = 1.83e-5;
den = 1.189;    
f = 1600;
Sample_Num = 28;
Desired_Freq = find (f == 1600);

% Defining X matrix
x = zeros(28,4);
x(1:7, 1) = [700e-6:50e-6:1000e-6];
x(8:28,1) = 700e-6;
x(1:7,2) = 1000e-6;
x(8:14, 2) = [700e-6:100e-6:1300e-6];
x(15:28,2) = 1000e-6;
x(1:14, 3) = 0.00615;
x(15:21, 3) = [0.004:0.001:0.010];
x(22:28,3) = 0.00615;
x(1:21, 4) = 0.0042;
x(22:28, 4) = [0.001:0.001:0.007];

y = zeros (Sample_Num, length(f));
for ii = 1 : 1 : Sample_Num
y(ii,:) = Sample_Calculation_3_edit (x(ii,:), ii, f, alpha, xi, c, n, den);
end
end

After that I used GA method in Optimization tool box to find the optimized x's (x(1), x(2), x(3) and x(4)) with boundary:

Lower boundary: [700e-6,700e-6,0.004,0.001]

Upper boundary: [1000e-6,1300e-6,0.010,0.007]

However, I got error "Your fitness function must return a scalar value.". The problem is i need to use all the value (scalar value) for optimization purpose.

What are the suitable optimization solver I can use to optimize the objective function with scalar value?

I tried to use multiobj, but it give me list of optimized parameter representing each sample available. But, I only want 1 set of optimized parameter (meaning that 1 optimized value for each x's).

Thank you.

Regards,

Nur Arafah

[Mario Malic]

I am not so good (if at all) at optimization stuff. I would suggest rather spend some time reading about optimization functions and think how your problem may relate to these. Maybe work with simpler problems from the optimization documentation and make your way up.

It kind of feels that with every message the problem expands in some sort of way or changes the objective, so try to be clear what do you want to do and explain it in the best way possible.

"I got error "Your fitness function must return a scalar value.". The problem is i need to use all the value (scalar value) for optimization purpose.", for GA, by reffering to documentation you can see that objective function must return a scalar value. About the other optimization algorithms that accept vectors see here LINK.

Also, I would recommend you to write the code for optimization rather than use it through toolbox. You'll understand more while doing so.

Here's a sample code.

ii = 1;
options = optimoptions('ga','TolFun', 1e-3, 'PlotFcn', @gaplotbestf);
lb = [700e-6,700e-6,0.004,0.001];
ub = [1000e-6,1300e-6,0.010,0.007];
[X_Opt, fval, exitflag,output] = ga(@(x)Sample_Calculation(x, ii, f, alpha, xi, c, n, den), 4,[], [], [], [], lb, ub, [], options)

And if you change the output in function to return the scalar value like this or in any other way the code will work.

yn = 1 - abs(RF).^2;  %% output for sample #1
y = 1 + yn(1,705); % obj fun

Sorry that I can not proceed further with the help, since I am not able to understand what you actually want to do.

[Jafar Nur Arafah]

Hi Mario Malic,

Thank you so much. You already give a lot inputs to me. I really appreciate it.

I will proceed accordingly to your advice.

Thank you again. Have a nice day.

Regards,

Nur Arafah

[Mario Malic]

You're welcome. Good luck and have a great day as well!