Fitting a model of a membrane with fminsearch

Hello guys,

I need to fit a model of a membrane with experimental datas and for that I am using fminsearch to calculate least squares. I have 15 inputs to fit.

The problem is I get this error:

my code for fitting is:

% Startwert
Lp = 10.9e-7;
%a=[7/60 -2227/50.400 4.0180 -3.9788 -1.9215 4.392 5.006];
a1=7/60;
a2=-2227/50.400;
a3=4.0180;
a4=-3.9788;
a5=-1.9215;
a6=4.392;
a7=5.006;
%n_a=numel(a)
%b=[29/20 97.453/50.400 -62.2680 -18.509 -12.156 9.15 43.17];
b1=29/20;
b2=97.453/50.400;
b3=-62.2680;
b4=-18.509;
b5=-12.156;
b6=9.15;
b7=43.17;
%n_b=numel(b)
%a_fit=zeros(1,n_a);
%b_fit=zeros(1,n_b);
initial=[a1 a2 a3 a4 a5 a6 a7 b1 b2 b3 b4 b5 b6 b7 Lp];

%Lp_fit=fminsearch(@error_function,Lp)
%fun=@(a,b,Lp)error_function(a,b,Lp);
a_fit=fminsearch(@error_function,initial)
%b_fit=fminsearch(@error_function,b)

and the error is:

Not enough input arguments.

Error in error_function (line 47)
Kt=(9/4)*(pi^2)*sqrt(2)*((1-lambda)^(-5/2))*(1 + ((a1)*(1-lambda)) - ((a2)*((1-lambda)^2))) + a3 - a4*lambda - (a5*(lambda)^2) +
(a6*(lambda)^3) + (a7*(lambda)^4);

Error in fminsearch (line 200)
fv(:,1) = funfcn(x,varargin{:});

Error in test (line 27)
a_fit=fminsearch(@error_function,initial)
 
>> 

can someone guide me?

Answer 1

[Matt J]

We can't see your "error_function", but you need to check its syntax and make sure it only expects a single input argument.

[shadan Kazempour Sabahi]

this is my error_function:

function error_function = f( a,b,Lp )
%%calculating Kc and Kd

Lp = 10.9e-7;                   % membrane hyraulic permeability [m/(s*bar)] 
a=[(7/60) -(2227/50.400) 4.0180 -3.9788 -1.9215 4.392 5.006];
%n=numel(a)
b=[(29/20) (97.453/50.400) -62.2680 -18.509 -12.156 9.15 43.17];

%Radius of the salute
r_s=1.83e-10; %[m]
%size the Membrane's pore
r_p_Nanomax=0.51e-9; %[m]
lambda=r_s/r_p_Nanomax;
%Kt
Kt=(9/4)*(pi^2)*sqrt(2)*((1-lambda)^(-5/2))*(1 + ((a(1))*(1-lambda)) - ((a(2))*((1-lambda)^2))) + a(3) - a(4)*lambda - (a(5)*(lambda)^2) + (a(6)*(lambda)^3) + (a(7)*(lambda)^4);

%Ks
Ks=(9/4)*(pi^2)*sqrt(2)*((1-lambda)^(-5/2))*(1 + ((b(1))*(1-lambda)) - ((b(2))*((1-lambda)^2))) - b(3) - b(4)*lambda - (b(5)*(lambda)^2) + (b(6)*(lambda)^3) + (b(7)*(lambda)^4);

%Kc and Kd
Kd = (6*pi)/Kt;
Kc = (1 + 2*lambda - lambda^2)*((Ks)/(2*Kt));

%parameters
V_hat_Feed = 4 / (1000 * 60);              % feed flow rate [m³/s]
%v_cs = 4 / (1000 * 60);         % crossflow rate [m/s]
A_R = 0.001; %V_hat_Feed / v_cs        % retentate channel cross-section assumed from feed flow and recommend v_cs [m²]
A_P = 10*A_R;                     % permeate channel cross-section [m²]  - assumption that it is a 10th of the retentate channel

A_M_total = 10*0.37;                        % membrane area [m²]
L_M = 1;                        % membrane module length inlet to outlet [m]
d_F1 = 0.001;                    % length of film on retentate side of membrane [m]
d_F2 = 0.001;                    % length of film on permeate side of membrane [m]
d_M = 150e-6;                   % membrane thickness [m]

% Salt
c0 = [4.3497 43.497 434.97];    % feed concentration [mol/m³]
%c0 = [4.3497 21.75 43.497 217.5 434.97];    % feed concentration [mol/m³]
num_c0=numel(c0);
D_Na = 1.33e-9;                 % diffusion coefficient of Na+ [m2/s]
D0_Na = 1.28e-9;                  % Diffusion coefficient of Na+ in very dilute solution  [m2/s]
M_Na = 22.98977;                  % molecular weight g/mol
r_Na = 0.98e-10;                  % size of Na
Tc = 10;                          % C
Tk = 273.15 + Tc;                 % K

%Factors for membrane
epsilon = 0.3;                  % porosity of the membrane
L_pore = 15e-7;                % length of the pore
                      
r_p=5e-10;                      % size of the pores [m]
phi=(1-(r_Na/r_p))^2;           % membrane porosity [-]

%Pressure
min_p=5;                % [bar]
max_p=40;               % [bar]
dp=1; 
p=min_p:dp:max_p;
num_p=numel(p);


%local discretization retentate & permeate (y)
y_min = 0;                  % [m]
y_max = L_M;                % [m]
y_steps = 5;                % process divided into x lumped capacites
dy = y_max / y_steps;       % [m]     -> value at y_max is only outlet conditions
yVec = y_min:dy:y_max;      
n_y = numel(yVec);

A_M = A_M_total / y_steps;
tau_y = A_R * dy / V_hat_Feed;
%local discretization film 1 (x)
x_min_F1 = 0;                  % [m]
x_max_F1 = d_F1;                  % [m]
x_steps_F1 = 20;                % number of discretization steps
dx_F1 = x_max_F1 / (x_steps_F1); % m     -> (x_steps) since 0 is only boundary conditions
xVec_F1 = x_min_F1:dx_F1:x_max_F1;
nx_F1 = numel(xVec_F1);
%local discretization film 2 (x)
x_min_F2 = 0;                  % [m]
x_max_F2 = d_F2;                  % [m]
x_steps_F2 = 20;                % number of discretization steps
dx_F2 = x_max_F2 / (x_steps_F2); % m     -> (x_steps) since 0 is only boundary conditions
xVec_F2 = x_min_F2:dx_F2:x_max_F2;
nx_F2 = numel(xVec_F2);
%local discretization membrane (x)
x_min_M = 0;                  % [m]
x_max_M = d_M;                % [m]
x_steps_M = 2;                % number of discretization steps
dx_M = x_max_M / (x_steps_M - 1); % m     -> (x_steps - 1) since 0 is first capacity
xVec_M = x_min_M:dx_M:x_max_M;
nx_M = numel(xVec_M);
%time discretization
dt = 0.5;                 % [s]
t_min = 0;              % [s]
t_max = 400;            % [s]
tVec = t_min:dt:(t_max+dt);
n_t = numel(tVec);



%vector set-up
c_i_R = zeros(n_y,n_t);
c_i_dot_R = zeros(n_y,n_t);
V_hat_R = zeros(n_y,n_t);

c_i_P = zeros(n_y,n_t);
c_i_dot_P = zeros(n_y,n_t);
V_hat_P = zeros(n_y,n_t);

c_i_F1 = zeros(n_y,nx_F1,n_t);
c_i_dot_F1 = zeros(n_y,nx_F1,n_t);

c_i_M = zeros(n_y,nx_M,n_t);
c_i_dot_M = zeros(n_y,n_t);

c_i_F2 = zeros(n_y,nx_F2,n_t);
c_i_dot_F2 = zeros(n_y,nx_F2,n_t);
n_i_hat_F2_out = zeros(n_y,n_t);

j_NaCl = zeros(n_y,n_t);
j = zeros(num_p,1);
R_obs = zeros(n_t,1);
R_obs_final = zeros(num_c0,num_p);


k_1_R = zeros(n_y,n_t);
k_2_R = zeros(n_y,n_t);
k_3_R = zeros(n_y,n_t);
k_4_R = zeros(n_y,n_t);

k_1_P = zeros(n_y,n_t);
k_2_P = zeros(n_y,n_t);
k_3_P = zeros(n_y,n_t);
k_4_P = zeros(n_y,n_t);

n_i_dot_total = zeros(n_t,1);
%initial conditions
c_i_R(:,1) = 0;                   % pure wate in the system at simulation start
c_i_F1(:,:,1) = 0;                % pure wate in the system at simulation start
c_i_M(:,:,1) = 0;                % pure wate in the system at simulation start
c_i_F2(:,:,1) = 0;                % pure wate in the system at simulation start
c_i_P(:,1) = 0;                   % pure wate in the system at simulation start
V_hat_P(n_y,:) = 0;
c_i_P(n_y,:) = 0;



%%SIMULATION SCRIPT

%Set-up
for i = 1:num_c0
    c_i_Feed = c0(i);
    
    for P = 1:num_p
        %S(i,P).cF1=[];
        %S(i,P).cF2=[];
        S(i,P).cR=[];
        S(i,P).cM=[];
        S(i,P).cP=[];
        S(i,P).R=[];
        S(i,P).ciRdot=[];
        S(i,P).ciPdot=[];
        S(i,P).nidot=[];
        
        j(P)=Lp*p(P); %Flow goes into the membrane
        tau_F1 = dx_F1 / j(P);
        tau_F2 = dx_F2 / j(P);

        for t = 1:n_t-1
%retentate R
% n_i_dot_R = n_i_hat_in - n_i_hat_F1 - n_i_hat_out
% n_i_dot_R = c_i_in * V_hat_in - c_i_R * j(p) * A_M - c_i_R * ( V_hat_in - j(p) * A_M )
% ...
% c_i_dot_R = V_hat_in * ( c_i_in - c_i_R ) / ( A_R * dy )
            for y = 1
                k_1_R(y,t) = V_hat_Feed * (c_i_Feed - c_i_R(y,t)) / (A_R * dy );        % Runge-Kutta Approach 4th Order (RK4)
                k_2_R(y,t) = V_hat_Feed * (c_i_Feed - (c_i_R(y,t) + (dt/2) * k_1_R(y,t))) / (A_R * dy );
                k_3_R(y,t) = V_hat_Feed * (c_i_Feed - (c_i_R(y,t) + (dt/2) * k_2_R(y,t))) / (A_R * dy );
                k_4_R(y,t) = V_hat_Feed * (c_i_Feed - (c_i_R(y,t) + dt * k_3_R(y,t))) / (A_R * dy );
                c_i_R(y,t+1) = c_i_R(y,t) + dt * ( (k_1_R(y,t)/6) + (k_2_R(y,t)/3) + (k_3_R(y,t)/3) + (k_4_R(y,t)/6) );
                
                V_hat_R(y,t) = V_hat_Feed - j(P) * A_M;
            end
            
            for y = 2:n_y-1
                k_1_R(y,t) = V_hat_Feed * (c_i_R(y-1,t) - c_i_R(y,t)) / (A_R * dy );
                k_2_R(y,t) = V_hat_Feed * ((c_i_R(y-1,t) + (dt/2) * k_1_R(y-1,t)) - (c_i_R(y,t) + (dt/2) * k_1_R(y,t))) / (A_R * dy );
                k_3_R(y,t) = V_hat_Feed * ((c_i_R(y-1,t) + (dt/2) * k_2_R(y-1,t)) - (c_i_R(y,t) + (dt/2) * k_2_R(y,t))) / (A_R * dy );
                k_4_R(y,t) = V_hat_Feed * ((c_i_R(y-1,t) + dt * k_3_R(y-1,t)) - (c_i_R(y,t) + dt * k_3_R(y,t))) / (A_R * dy );
                c_i_R(y,t+1) = c_i_R(y,t) + dt * ( (k_1_R(y,t)/6) + (k_2_R(y,t)/3) + (k_3_R(y,t)/3) + (k_4_R(y,t)/6) );
                
                V_hat_R(y,t) = V_hat_R(y-1,t) - j(P) * A_M;
            end
            
            for y = n_y                         % since last step is only outflow; y+dy is one capacity
                c_i_R(y,t+1) = c_i_R(y-1,t+1);
                
                V_hat_R(y,t) = V_hat_R(y-1,t);
            end
            
            c_i_dot_R(:,t) = ( c_i_R(:,t+1) - c_i_R(:,t) ) / dt;        % estimated just for plotting
            
%film F1
% n_i_dot_F1 = n_i_hat_in - n_i_hat_out
% c_i_dot_F1 * V_F1 = c_i_in * j(p) * A_M + A_M * ( D_Na * dc_i/dx - c_i_out * v_x )
% ... dc_i/dx = (c_i_out - c_i_in)/dx_F1 ...
% c_i_dot_F1 = c_i_in * j(p) / dx_F1  +  (D_Na / dx_F1) * (c_i_out - c_i_in) / dx_F1  -  c_i_out * j(p) / dx_F1

%%% KEEP IN MIND: y = n_y no film diffusion since @ n_y is just outlet of the membrane module
%             c_i_F1(:,1,t) = c_i_R(:,t);
%             c_i_dot_F1(:,1,t) = 0;
%             
%             for x_F1 = 1:nx_F1-2
%                 c_i_F1(:,x_F1+1,t+1) = c_i_F1(:,x_F1+1,t) + dt*((D_Na/((dx_F1)^2))*(c_i_F1(:,x_F1+2,t)-2*c_i_F1(:,x_F1+1,t)+c_i_F1(:,x_F1,t)) - j(P) * (c_i_F1(:,x_F1+2,t) - c_i_F1(:,x_F1,t))/(2*dx_F1));  %FDM
%             end
%             
%             for x_F1 = nx_F1
%                 c_i_F1(:,x_F1,t+1) = c_i_F1(:,x_F1,t) + dt*((D_Na/((dx_F1)^2))*(c_i_F1(:,x_F1,t)-2*c_i_F1(:,x_F1-1,t)+c_i_F1(:,x_F1-2,t)) - j(P) * (c_i_F1(:,x_F1,t)-c_i_F1(:,x_F1-1,t)) / (dx_F1));      %Backward FDM
%             end
%              
%membrane M
            c_i_M(:,1,t) = phi * c_i_R(:,t);
            c_i_M(:,nx_M,t) = phi * c_i_P(:,t);
            Pe_pore=(Kc*j(P)*L_pore)/(Kd*D0_Na*epsilon); % Peclet Number in the pores
            
            for y = 1:n_y-1
                j_NaCl(y,t) = (Kc*j(P)/epsilon) * ((c_i_M(y,nx_M,t) - c_i_M(y,1,t)*(exp(Pe_pore)))/(1-exp(Pe_pore)));
            end
  
%permeate
 
for y = n_y-1:-1:1              % backwards loop since conditions in P @ y+1 have to be known before estimating P @ y  -  n_y would only be not zero if a feed to P exists  -  k-values if n_y are zero anyway
                k_1_P(y,t) = ( c_i_P(y+1,t) * V_hat_P(y+1,t) +  j_NaCl(y,t) * A_M - c_i_P(y,t) * ( V_hat_P(y+1,t) + j(P) * A_M )  ) / (A_P * dy);
                k_2_P(y,t) = ( (c_i_P(y+1,t) + (dt/2) * k_1_P(y+1,t)) * V_hat_P(y+1,t) +  j_NaCl(y,t) * A_M - (c_i_P(y,t) + (dt/2) * k_1_P(y,t)) * ( V_hat_P(y+1,t) + j(P) * A_M )  ) / (A_P * dy);
                k_3_P(y,t) = ( (c_i_P(y+1,t) + (dt/2) * k_2_P(y+1,t)) * V_hat_P(y+1,t) +  j_NaCl(y,t) * A_M - (c_i_P(y,t) + (dt/2) * k_2_P(y,t)) * ( V_hat_P(y+1,t) + j(P) * A_M )  ) / (A_P * dy);
                k_4_P(y,t) = ( (c_i_P(y+1,t) + (dt) * k_4_P(y+1,t)) * V_hat_P(y+1,t) +  j_NaCl(y,t) * A_M - (c_i_P(y,t) + (dt) * k_4_P(y,t)) * ( V_hat_P(y+1,t) + j(P) * A_M )  ) / (A_P * dy);
                c_i_P(y,t+1) = c_i_P(y,t) + dt * ( (k_1_P(y,t)/6) + (k_2_P(y,t)/3) + (k_3_P(y,t)/3) + (k_4_P(y,t)/6) );
                
                V_hat_P(y,t) = V_hat_P(y+1,t) + j(P) * A_M;
                c_i_dot_P(y,t) = ( c_i_P(y,t+1) - c_i_P(y,t) ) / dt;        % estimated just for plotting
            end


%Finish-up
            n_i_dot_total(i,t) = c0(i) * V_hat_Feed - c_i_P(1,t) * V_hat_P(1,t) - c_i_R(n_y,t) * V_hat_R(n_y,t);
            R_obs(t) = 1-(c_i_P(1,t)/c0(i));
            
        end 

        R_obs_final(i,P) = 1-(c_i_P(1,n_t-1)/c0(i));
        
        %S(i,P).cF1=c_i_F1;
        %S(i,P).cF2=c_i_F2;
        S(i,P).cR=c_i_R;
        S(i,P).cM=c_i_M;
        S(i,P).cP=c_i_P;
        S(i,P).R=R_obs;
        S(i,P).ciRdot=c_i_dot_R;
        S(i,P).ciPdot=c_i_dot_P;
        S(i,P).nidot=n_i_dot_total;

    end 
end

%Real Data
L_permeate=[438 350 337 335 332 330];
n=numel(L_permeate);
L_Feed=[1281 1279 1278 1280 1281 1280];
n=numel(L_Feed);
L_retentate=[1294 1320 1331 1340 1351 1375];
p2=[10 15 20 25 30 40];
n_p=numel(p2);
R_obs1=zeros(n,1);

for i=1:n
    R_obs1(i)=1-(L_permeate(i)/L_Feed(i));
   
end

%%Show the results
figure(1)
plot(p2,R_obs1,'r*-',p(5:35),R_obs_final(1,5:35))
legend('real data','simulation results','Location','SouthEast')
grid on
title('Retention vs. Transmembrane pressure')
xlabel('transmembrane pressure [bar]')
ylabel('R_o_b_s')

%fitting coefficients
% R_obs1=R_obs1';
% polyfit(p2,R_obs1,3);
R=[R_obs_final(1,5) R_obs_final(1,10) R_obs_final(1,15) R_obs_final(1,20) R_obs_final(1,25) R_obs_final(1,35)];
R=R';

error_function=sum((R-R_obs1).^2);

end

[Matt J]

You need to wrap it in a single-argument function.

fun=@(x) error_function(x(1:7), x(8:14),x(15));

a_fit=fminsearch(fun,initial)

[shadan Kazempour Sabahi]

thank you for the answer.

[shadan Kazempour Sabahi]

Now a new problem apeared. fminsearch always gives me the same start value. not the optimized one, exactly the same value that I give as the start value. What could be the reason for that?