Optimization with definite integration

I have a cost function (sai) which I need to optimize to minimize the difference between three signals (namely power spectral densities).

Syy, S1y and S2y are all same size and are powers spectral and cross spectral densities respectively. L is an integer and stands for length range of integration (ohm) is from zero to 2-3 (depending the on the first peak of frequency) H1_bar and H2bar are as follows:

here is how i defined my optimization function:

if true
    % error_norm = @ (x) integral (@ (sai)(Syy -  (((0.5*(((sin (sai)+ sinh(sai))/((sin(sai)*cosh(sai))+(cos(sai)*sinh(sai)))))) .* S1y) + ...
(((-1 / (4 * sai)) *(((cos (sai)- cosh(sai))/((sin(sai)*cosh(sai))+(cos(sai)*sinh(sai))))) * sub_length) .* S2y))) .^ 2,0,x) ;

%[x,fval,exitflag,output] = fmincon (error_norm,x0,A,b,Aeq,beq,lb,ub,nonlcon);
  end

I am getting error that inner matrix dimension must agreee. When I delete the integration I receive this error (user input value must return a scalar).

1. I really deeply appreciate if someone help me out, I am stuck in this part for almost a week
2. Also, how may I plot this so I can have an idea what would be the approximate initial guess for optimization.

Thanks ahead :)

Answer 1

[Torsten]

1. Integrate the following functions over big_Omega using MATLAB's "trapz":

Syy*S1y
Syy*S2y
(S1y)^2
S1y*S2y
(S2y)^2.

Let the results be

I1
I2
I3
I4
I5

2. Compute dh1_bar/dsai, dh2_bar/dsai (symbolically).

3. Solve the equation

I1*dh1_bar/dsai + I2*L*dh2_bar/dsai - I3*h1_bar*dh1_bar/dsai - I5*L^2*h2_bar*dh2_bar/dsai - I4*L*(h1_bar*dh2_bar/dsai + h2_bar*dh1_bar/dsai)) = 0

for sai using MATLAB's "fzero".

The sai you get is optimal.

The method used is to simply differentiate "under the integral" with respect to sai and to set the derivative to zero.

Best wishes

Torsten.

[Shz713]

Dear Torsten,

Heaps of thanks for your insight. I'm away from my PC but I will do as you instructed. Is it okay I post my results here tommorrow for you to check? I'm grateful if you could follow-up my post.

Cheers

[Shz713]

Dear Torsten, check out the code below: 1- How may I integrate over big omega using trapz? should I define line spaicng or just integrate over the points of interest?

2- when I run the fzero I get this error "Function values at interval/initial guess endpoints must be finite and real." What should be ze the x0 value?

if true
  % D_1 = load ('d_1.txt');                 %vertical acceleration @ initial
D_3 = load ('d_3.txt');                 %vertical acceleration @ end
Fi_1 = load ('fi_1.txt');               %rotational acceleration @ initial
Fi_3 = load ('fi_3.txt');               %rotational acceleration @ end
D_2 = load ('d2.txt');                  % vertical acceleration at middle
u_1 = D_1 + D_3;
u_2 = Fi_1 - Fi_3;             
L = 39.3;   
YY = abs (fft (D_2));                    %FFT of signal
Y_1 = abs (fft (u_1));                   %abs gives magnitude spectrum, otherwise shows compelx numbers
Y_2 = abs (fft (u_2));                        %by default bins are added by one in Matlab->(N-1)bins are started from zero not one
Syy= YY.* conj(YY);                       % power spectral density (PSD)  
S1y = Y_1.*conj (YY);                     % cross spectral density (CSD)
S2y = Y_2.*conj (YY);
spacing = 0:0.1:100;
G1 = Syy .* S1y;
I1 = trapz(G1);      %I1 = trapz(G1(1:100));
G2 = Syy .* S2y;
I2 = trapz (G2);
G3 = S1y .^2;
I3 = trapz (G3);
G4 = S1y .* S2y;
I4 = trapz(G4);
G5 = S2y .^2;
I5 = trapz (G5);
syms sai
h1_bar = 0.5*((sin (sai)+ sinh(sai))/((sin(sai)*cosh(sai))+(cos(sai)*sinh(sai))));
h2_bar = ((-1 / (4 * sai)) *((cos (sai)- cosh(sai))/((sin(sai)*cosh(sai))+(cos(sai)*sinh(sai)))));
dh1_bar = diff (h1_bar);
dh2_bar = diff (h2_bar);
mysai = @ (sai) (I1 * dh1_bar) + (I2 * L * dh2_bar) - (I3 * h1_bar * dh1_bar) - (I5 * L^2 * h2_bar * dh2_bar) -...
        (I4 * L *((h1_bar * dh2_bar) + (h2_bar * dh1_bar)));

endx0=[0 9];
zeta = fzero (mysai,x0);

Thanks again for your assistance

[Torsten]

1. When calculating the integrals, you will have to supply the corresponding omega values, too. E.g. I1=trapz(omega,G1).

2. Since fzero can not handle symbolic expressions, you will have to transform dh1_bar and dh2_bar into a function handle. Also define h1_bar and h2_bar as function handles.

3. First try to plot mysai and see where the zero is located before calling fzero.

Best wishes

Torsten.

Answer 2

[Torsten]

Why do you integrate with respect to big_omega whereas in the integral you use small_omega ?

What is big_omega ? An interval ?

How do you supply S_yy, S_1y, S_2y ? As a table of values ? As an explicit function ?

sai is a single scalar value ? h1_bar(sai), h2_bar(sai) are single scalar values ?

Best wishes

Torsten.

[Shz713]

Thanks for your prompt response mate,

I did not understand your first Q, omega is to indicate the parameters must be in frequency domain not time domain, hence integration is with respect to ohm, which is the range of first natural frequency based

2-Syy, S1y and S2y are powers spectral densities, or in other words time signals converted to frequency domain by pwelch/fft; hence it gives column vector depending on size of data. In my case is [501,1]. Up until the minimization, these functions are running smoothly with correct graph. So all of them have same size

3-indeed, Sai must be scalar value to minimize the difference of Syy, S1y and S2y. Similarly, h1_bar and h2_bar will be scalar values.

The main concept is to minimize the difference between Syy- (S1y+S2y). So h1_bar and h2_bar are multiplied as objective parameters to minimize this difference.

Do you want me to post the entire code to be easier to catch up? In fact this is algorithm defined in a journal article n I can attach that as well.

Warm regards Sir