Why optimization has a Initial point value

Why optimization has a Initial point value

[Geoff Hayes]

Pearwpun - I think most optimization algorithms have an initial point so that the algorithm can then improve upon that point.

Why do you think that it is unnecessary?

[Pearwpun Bunjun]

I means answer is from my code

[Geoff Hayes]

Oh you mean why is the answer the same as your initial value? Is that your question?

Why did you choose 310 as an initial point?

[Pearwpun Bunjun]

it's initial of process

Answer 1

[Pearwpun Bunjun]

    X0    = 310*ones(1,m);
    Lb   =(30+273)*ones(1,m);
    Ub    =(50+273)*ones(1,m);
    nonlcon = @constraints;
   options = optimoptions(@fmincon,'Algorithm','sqp','MaxIter',300);
   [MV,fval]= fmincon(@Obj,X0,[],[],[],[],[],[],@NONLCON,options);

answer is X0 = 310

[Torsten]

What is your question ?

[Pearwpun Bunjun]

why answer from fmincon has X0 Single value

[Torsten]

What is @Obj ?

What is @NONLCON (note that "nonlcon" and "NONLCON" are different because MATLAB is case-sensitive) ?

[Pearwpun Bunjun]

@Obj=fun @NONLCON=constraints correct @nonlcon is model incorrect did not use

[Torsten]

We have to see these two functions in order to answer your question.

[Pearwpun Bunjun]

function PPP7 clear global C1 Tj1 Tjsp1 p m

b = 1.45; % dimensionless : nucleation rate exponential g = 1.5; % dimensionless : growth rate exponential kb = 285.0; % the nucleation rate constant (1/(s um3) kg = 1.44*10^8; % the growth rate constant (um/s) EbR = 7517.0; % Eb/R : the nucleation activation energy/gas constant (K) EgR = 4859.0; % Eg/R : the growth activation energy/gas constant(K) U = 1800; % the overall heat transfer coefficient (kJ/(m2 h K)) A = 0.25; % the total heat transfer surface area (m2) delH = 44.5; % the heat of reaction (kJ/kg) Cp = 3.8; % the heat capacity of the solution (kJ/(kg K)) M = 27.0; % the mass of solvent in the crystallizer (kg) rho = 2.66*10^-12; % the density of crystal (g/um3) kv = 1.5; % the volumetric shape factor

%jacket Parameter Vj=0.015; %m3 Fj=0.001; %m3/s rhoj=1000; %kg/m3 Cpj=4.184; %J/kgK

% Time TTime=30; tf = (60*TTime); % sec (30 min) dt = 20; % sec nt = tf/dt; t = linspace(0,TTime,nt+1);% min sampt= 2; % unit: min % q=1; p=30; m=30; z=0; % initial conditions

C(1) = 0.1743; % g solute/g solvent T(1) = 323; % K Tj(1) =278.3; % K Tjsp(1) = (20+273); % K initial values of manipulated variable %Tsp(1) =(127+273);

%C1(1)=C(1); %T1(1)=T(1); %Tj1(1)=Tj(1); %Tsp1(1)=Tsp(1);

%for Y=1:90

%Tsp(20)=490; % Tsp(Y+1)=Tsp(Y);

%Tsp(52)=450; % Tsp(Y+1)=Tsp(Y); %Tsp(70)=430; % Tsp(Y+1)=Tsp(Y); %end

% Parabolic distribution @ 70-90 um un0(1) = 0; un1(1) = 0; un2(1) = 0; un3(1) = 0; un4(1) = 0; un5(1) = 0;

us0(1) = 70;

us1(1) = 1.8326e+004; us2(1) = 5.0480e+006; us3(1) = 1.3928e+009; us4(1) = 3.8490e+011; us5(1) = 1.0654e+014;

h=waitbar(0,'Simulation in Process...');

for i=1:nt; waitbar(i/nt); % Moment models % the total crystal number uu0(i) = un0(i)+us0(i); % the total crystal length uu1(i) = un1(i)+us1(i); % the total crystal surface area uu2(i) = un2(i)+us2(i); % the total crystal volume uu3(i) = un3(i)+us3(i);

    % Saturation concentration (T in celcious)
    Temp(i) = T(i)-273;
    Cs(i) = 6.29*10^-2+(2.46*10^-3*Temp(i))-(7.14*10^-6*Temp(i)^2);
    % Metastable concentration
    Cm(i) = 7.76*10^-2+(2.46*10^-3*Temp(i))-(8.10*10^-6*Temp(i)^2);

    % Supersaturation
    S(i) = (C(i)-Cs(i))/Cs(i);
    % The nucleation rate
    B(i) = kb*(exp(-EbR/T(i)))*(((C(i)-Cs(i))/Cs(i))^b)*uu3(i);

    % The growth rate
    G(i) = kg*(exp(-EgR/T(i)))*(((C(i)-Cs(i))/Cs(i))^g);

    % Population Balance Equation(PBE) 

     nr=600;
     for j=1:nr+1  % r=0:600
        r = j+1;

     % Seed
     % Initial condition
       if r>=250 & r<=300 
          rs(j,1) = r; 

     % Parabolic distribution
          ns(j,1) = 0.0032*(300-r)*(r-250);
       else
          rs(j,1) = r;
          ns(j,1) = 0;
       end  %  if j>=250 &j<=300 
    rs(j,i+1) = rs(j,i)+G(i)*dt;
    ns(j,i+1) = ns(j,i);  

    if rs(j,i+1)>600
        rs(j,i+1)=600;
    end

      % Nucleation
        % at t=0
       if i==1
            rn(1,1) = 0;
            nn(1,1) = B(1)/G(1);

       else % at t=dt
           if j==1 % r=0
               rn(1,i) = 0;
               nn(1,i) = B(i)/G(i);
           else
               rn(j,i) = rn(j-1,i-1)+G(i)*dt;
               nn(j,i) = nn(j-1,i-1);  
           end % if j==1
       end % if i==1      
        end % for j=1:nr+1

   % The crystal size distribution  

    %%%%%%%%The crystal size distribution
    if i==1
        Nn = [B(1)/G(1); zeros(nr,1)];
    else
        Nn = nn([1:nr+1],i);
    end
        Ns = ns([1:nr+1],i);

        n = Nn+Ns;

    if i==1
        rn([1:nr+1],i) = 0;
    else
       rn(nr+1,i) = rn(nr,i);
    end
       rrn = rn([1:nr+1],i);
       rrs = rs([1:nr+1],i);

   % Moment : Trapezoidal Rule

   for k=1:nr+1

     % u2
        if k==1
        sumun22 = rrn(1,1)^2*Nn(1,1);
        sumus22 = rrs(1,1)^2*Ns(1,1);
        else
        sumun22 = sumun22+(((rrn(k-1,1)^2*Nn(k-1,1))+(rrn(k,1)^2*Nn(k,1)))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus22 = sumus22+(((rrs(k-1,1)^2*Ns(k-1,1))+(rrs(k,1)^2*Ns(k,1)))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
        sumu2=sumun22+sumus22;

    % u3
        if k==1
        sumun33 = rrn(1,1)^3*Nn(1,1);
        sumus33 = rrs(1,1)^3*Ns(1,1);
        else
        sumun33 = sumun33+(((rrn(k-1,1)^3*Nn(k-1,1))+(rrn(k,1)^3*Nn(k,1)))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus33 = sumus33+(((rrs(k-1,1)^3*Ns(k-1,1))+(rrs(k,1)^3*Ns(k,1)))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
            sumu3=sumun33+sumus33;

    % u1
        if k==1
        sumun11 = rrn(1,1)^1*Nn(1,1);
        sumus11 = rrs(1,1)^1*Ns(1,1);
        else
        sumun11 = sumun11+(((rrn(k-1,1)^1*Nn(k-1,1))+(rrn(k,1)^1*Nn(k,1)))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus11 = sumus11+(((rrs(k-1,1)^1*Ns(k-1,1))+(rrs(k,1)^1*Ns(k,1)))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
            sumu1=sumun11+sumus11; 

     % u0
        if k==1
        sumun00 = rrn(1,1)^0*Nn(1,1);
        sumus00 = rrs(1,1)^0*Ns(1,1);
        else
        sumun00 = sumun00+((rrn(k-1,1)^0*Nn(k-1,1)+rrn(k,1)^0*Nn(k,1))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus00 = sumus00+((rrs(k-1,1)^0*Ns(k-1,1)+rrs(k,1)^0*Ns(k,1))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
            sumu0=sumun00+sumus00;          

     % u4
        if k==1
        sumun44 = rrn(1,1)^4*Nn(1,1);
        sumus44 = rrs(1,1)^4*Ns(1,1);
        else
        sumun44 = sumun44+(((rrn(k-1,1)^4*Nn(k-1,1))+(rrn(k,1)^4*Nn(k,1)))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus44 = sumus44+(((rrs(k-1,1)^4*Ns(k-1,1))+(rrs(k,1)^4*Ns(k,1)))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
            sumu4=sumun44+sumus44;

        % u5
        if k==1
        sumun55 = rrn(1,1)^5*Nn(1,1);
        sumus55 = rrs(1,1)^5*Ns(1,1);
        else
        sumun55 = sumun55+(((rrn(k-1,1)^5*Nn(k-1,1))+(rrn(k,1)^5*Nn(k,1)))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus55 = sumus55+(((rrs(k-1,1)^5*Ns(k-1,1))+(rrs(k,1)^5*Ns(k,1)))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
            sumu5=sumun55+sumus55;

    end % for k
    u2(i) = sumu2;
    u3(i) = sumu3;
    u1(i) = sumu1;
    u0(i) = sumu0;  
    u4(i) = sumu4;
    u5(i) = sumu5;

C1=C(i); Tjsp1=Tjsp(i); Tj1=Tj(i);

    % Mass Balance : Solute concentration
    C(i+1) = C(i)+dt*(-3*rho*kv*G(i)*u2(i));

    % Batch Energy Balance
    TT3=T(i);

    RR1=(-U*A/(M*Cp*3600)*(T(i)-Tj(i))-delH/Cp*3*rho*kv*G(i)*u2(i));
    T(i+1)=T(i)+dt*RR1;
    TT4=T(i+1);
    Temp(i+1)= T(i+1)-273;
    Cs(i+1) = 6.29*10^-2+(2.46*10^-3*Temp(i+1))-(7.14*10^-6*Temp(i+1)^2);
    Cm(i+1) = 7.76*10^-2+(2.46*10^-3*Temp(i+1))-(8.10*10^-6*Temp(i+1)^2);

    Tj(i+1)= Tj(i)+dt*(Fj/Vj*(Tjsp(i)-Tj(i))+(((U*A/3600)*(T(i)-Tj(i))/rhoj*Vj*Cpj)));
    Tjsp(i+1)=Tjsp(i);

    % Part III: Calculating the value of manipulated variable  % 

    X0    = T(1)*ones(1,m);
    Lb   =(30+273)*ones(1,m);
    Ub    =(50+273)*ones(1,m);
    nonlcon = @constraints;
   options = optimoptions(@fmincon,'Algorithm','sqp','MaxIter',300);
  % options = optimset('Display','iter');
   [MV,fval]= fmincon(@Obj,X0,[],[],[],[],[],[],@NONLCON,options);
   %x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)

    %@(MV)NONLCON(MV)
    SS4=MV;
    SS1=MV(1);
    T(i+1) = MV(1);
    %Moment model
    %u3
         %Neacleated class
            un3(i+1) = sumun33+dt*(3*G(i)*sumun22);
         %Seed class
         %Virtual process
            us3(i+1) = sumus33+dt*(3*G(i)*sumus22);
         %the total crystal volume
            uu3(i+1)=un3(i+1)+us3(i+1);
    %u2
         %Neacleated class
            un2(i+1) = sumun22+dt*(2*G(i)*sumun11);
         %Seed class
         %Virtual process
            us2(i+1) = sumus22+dt*(2*G(i)*sumus11);
         %the total crystal surface area
            uu2(i+1)=un2(i+1)+us2(i+1);

    %u1
         %Neacleated class
            un1(i+1) = sumun11+dt*(1*G(i)*sumun00);
         %Seed class
         %Virtual process
            us1(i+1) = sumus11+dt*(1*G(i)*sumus00);
         %the total crystal volume
            uu1(i+1)=un1(i+1)+us1(i+1);

    %u0
         %Neacleated class
            un0(i+1) = sumun00+dt*(0*G(i)*sumun00);
         %Seed class
         %Virtual process
            us0(i+1) = sumus00+dt*(0*G(i)*sumus00);
         %the total crystal volume
            uu0(i+1)=un0(i+1)+us0(i+1);

end delete(h) % DELETE the waitbar; don't try to CLOSE it.

save ppp7.mat RR1 TT3 TT4 figure subplot(2,1,1) stairs(t,T,'-b') legend('Reactor Temp') ylabel('T(k)') xlabel('Time (min)')

subplot(2,1,2) plot(t,C,'-b') legend('Concentration') xlabel('Time (min)') ylabel('C')

figure subplot(2,1,1) plot(t,Tj,'r',t,Tjsp,'b') legend('Tj') xlabel('Time (min)') ylabel('Tj')

figure subplot(1,1,1)

stairs(t,Tjsp,'r') legend('Tjsp') xlabel('Time (min)') ylabel('Tjsp')

function f = Obj(MV,m,p,C,Tj,Tjsp)

global C1 Tj1 Tjsp1 p m Lw II1 II2 II3 II4 II5

% Process parameters : Potassium sulfate (K2SO4-H2O) b = 1.45; % dimensionless : nucleation rate exponential g = 1.5; % dimensionless : growth rate exponential kb = 285.0; % the nucleation rate constant (1/(s um3) kg = 1.44*10^8; % the growth rate constant (um/s) EbR = 7517.0; % Eb/R : the nucleation activation energy/gas constant (K) EgR = 4859.0; % Eg/R : the growth activation energy/gas constant(K) U = 1800; % the overall heat transfer coefficient (kJ/(m2 h K)) A = 0.25; % the total heat transfer surface area (m2) delH = 44.5; % the heat of reaction (kJ/kg) Cp = 3.8; % the heat capacity of the solution (kJ/(kg K)) M = 27.0; % the mass of solvent in the crystallizer (kg) rho = 2.66*10^-12; % the density of crystal (g/um3) kv = 1.5; % the volumetric shape factor %jacket Parameter Vj=0.015; %m3 Fj=0.001; %m3/s rhoj=1000; %kg/m3 Cpj=4.184; %J/kgK

% Step size, Sampling time and process time % tf = 1800; % sec (30 min) dt = 20; % sec nt = tf/dt; t = linspace(0,30,nt+1);% min sampt= 2; % unit: min %

Cobj = C1 ; % g solute/g solvent Tobj = MV ; % K Tjobj = Tj1 ; % K Tjspobj = Tjsp1;

q=1; z=0;

% Parabolic distribution @ 70-90 um un0(1) = 0; un1(1) = 0; un2(1) = 0; un3(1) = 0; un4(1) = 0; un5(1) = 0;

us0(1) = 70;

us1(1) = 1.8326e+004; us2(1) = 5.0480e+006; us3(1) = 1.3928e+009; us4(1) = 3.8490e+011; us5(1) = 1.0654e+014;

for C=1:1:p

    if C>m
      Tobj(C)=  Tobj(C-1) 
    end
  end   

for  z=1:p

    % Moment models
    % the total crystal number
    uu0(z) = un0(z)+us0(z);
    % the total crystal length
    uu1(z) = un1(z)+us1(z);
    % the total crystal surface area
    uu2(z) = un2(z)+us2(z);
    % the total crystal volume
    uu3(z) = un3(z)+us3(z);

    % Saturation concentration (T in celcious)
    Temp(z) = Tobj(z)-273;
    Cs(z) = 6.29*10^-2+(2.46*10^-3*Temp(z))-(7.14*10^-6*Temp(z)^2);
    % Metastable concentration
    Cm(z) = 7.76*10^-2+(2.46*10^-3*Temp(z))-(8.10*10^-6*Temp(z)^2);

    % Supersaturation
    S(z) = (Cobj(z)-Cs(z))/Cs(z);
    % The nucleation rate
    B(z) = kb*(exp(-EbR/Tobj(z)))*(((Cobj(z)-Cs(z))/Cs(z))^b)*uu3(z);

    % The growth rate
    G(z) = kg*(exp(-EgR/Tobj(z)))*(((Cobj(z)-Cs(z))/Cs(z))^g);

     % Population Balance Equation(PBE) 

     nr=600;
     for j=1:nr+1  % r=0:600
        r = j+1;

     % Seed
     % Initial condition
       if r>=250 & r<=300 
          rs(j,1) = r; 

     % Parabolic distribution
          ns(j,1) = 0.0032*(300-r)*(r-250);
       else
          rs(j,1) = r;
          ns(j,1) = 0;
       end  %  if j>=250 &j<=300 
    rs(j,z+1) = rs(j,z)+G(z)*dt;
    ns(j,z+1) = ns(j,z);  

    if rs(j,z+1)>600
        rs(j,z+1)=600;
    end

      % Nucleation
        % at t=0
       if z==1
            rn(1,1) = 0;
            nn(1,1) = B(1)/G(1);

       else % at t=dt
           if j==1 % r=0
               rn(1,z) = 0;
               nn(1,z) = B(z)/G(z);
           else
               rn(j,z) = rn(j-1,z-1)+G(z)*dt;
               nn(j,z) = nn(j-1,z-1);  
           end % if j==1
       end % if i==1      
        end % for j=1:nr+1

   % The crystal size distribution  

    %%%%%%%%The crystal size distribution
    if z==1
        Nn = [B(1)/G(1); zeros(nr,1)];
    else
        Nn = nn([1:nr+1],z);
    end
        Ns = ns([1:nr+1],z);

        n = Nn+Ns;

    if z==1
       rn([1:nr+1],z) = 0;
    else
       rn(nr+1,z) = rn(nr,z);
    end
       rrn = rn([1:nr+1],z);
       rrs = rs([1:nr+1],z);

   % Moment : Trapezoidal Rule

   for k=1:nr+1

     % u2
        if k==1
        sumun22 = rrn(1,1)^2*Nn(1,1);
        sumus22 = rrs(1,1)^2*Ns(1,1);
        else
        sumun22 = sumun22+(((rrn(k-1,1)^2*Nn(k-1,1))+(rrn(k,1)^2*Nn(k,1)))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus22 = sumus22+(((rrs(k-1,1)^2*Ns(k-1,1))+(rrs(k,1)^2*Ns(k,1)))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
        sumu2=sumun22+sumus22;

    % u3
        if k==1
        sumun33 = rrn(1,1)^3*Nn(1,1);
        sumus33 = rrs(1,1)^3*Ns(1,1);
        else
        sumun33 = sumun33+(((rrn(k-1,1)^3*Nn(k-1,1))+(rrn(k,1)^3*Nn(k,1)))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus33 = sumus33+(((rrs(k-1,1)^3*Ns(k-1,1))+(rrs(k,1)^3*Ns(k,1)))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
            sumu3=sumun33+sumus33;

    % u1
        if k==1
        sumun11 = rrn(1,1)^1*Nn(1,1);
        sumus11 = rrs(1,1)^1*Ns(1,1);
        else
        sumun11 = sumun11+(((rrn(k-1,1)^1*Nn(k-1,1))+(rrn(k,1)^1*Nn(k,1)))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus11 = sumus11+(((rrs(k-1,1)^1*Ns(k-1,1))+(rrs(k,1)^1*Ns(k,1)))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
            sumu1=sumun11+sumus11; 

     % u0
        if k==1
        sumun00 = rrn(1,1)^0*Nn(1,1);
        sumus00 = rrs(1,1)^0*Ns(1,1);
        else
        sumun00 = sumun00+((rrn(k-1,1)^0*Nn(k-1,1)+rrn(k,1)^0*Nn(k,1))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus00 = sumus00+((rrs(k-1,1)^0*Ns(k-1,1)+rrs(k,1)^0*Ns(k,1))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
            sumu0=sumun00+sumus00;          

    % u4
        if k==1
        sumun44 = rrn(1,1)^3*Nn(1,1);
        sumus44 = rrs(1,1)^3*Ns(1,1);
        else
        sumun44 = sumun44+(((rrn(k-1,1)^3*Nn(k-1,1))+(rrn(k,1)^3*Nn(k,1)))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus44 = sumus44+(((rrs(k-1,1)^3*Ns(k-1,1))+(rrs(k,1)^3*Ns(k,1)))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
            sumu4=sumun44+sumus44;      

        % u5
        if k==1
        sumun55 = rrn(1,1)^3*Nn(1,1);
        sumus55 = rrs(1,1)^3*Ns(1,1);
        else
        sumun55 = sumun55+(((rrn(k-1,1)^3*Nn(k-1,1))+(rrn(k,1)^3*Nn(k,1)))*(rrn(k,1)-rrn(k-1,1))/2);
        sumus55 = sumus55+(((rrs(k-1,1)^3*Ns(k-1,1))+(rrs(k,1)^3*Ns(k,1)))*(rrs(k,1)-rrs(k-1,1))/2);
        end % if k==1
            sumu5=sumun55+sumus55;           

   end       

    u2(z) = sumu2;
    u3(z) = sumu3;
    u1(z) = sumu1;
    u0(z) = sumu0;  
    u4(z) = sumu4;
    u5(z) = sumu5;

% Mass Balance : Solute concentration
     TT1=Tobj(z);
    RR5 =real(-3*rho*kv*G(z)*u2(z));
    Cobj(z+1) = Cobj(z)+dt*RR5;
% Batch Energy Balance
    WW=real(-U*A/(M*Cp*3600)*(Tobj(z)-Tjobj(z))-delH/Cp*3*rho*kv*G(z)*u2(z));
    Tobj(z+1)= Tobj(z)+dt*WW;
    TT2=Tobj(z+1);
    Temp(z+1)= Tobj(z+1)-273;
    Cs(z+1) = (6*10^-5)*exp(0.0396*Temp(z+1));
    Tjobj(z+1)= Tjobj(z)+dt*(Fj/Vj*(Tjspobj(z)-Tjobj(z))+(((U*A/3600)*(Tobj(z)-Tjobj(z))/rhoj*Vj*Cpj)));
    Tjspobj(z+1)=Tjspobj(z);

   %Moment model
    %u3
         %Neacleated class
            un3(z+1) = sumun33+dt*(3*G(z)*sumun22);
         %Seed class
         %Virtual process
            us3(z+1) = sumus33+dt*(3*G(z)*sumus22);
         %the total crystal volume
            uu3(z+1)=un3(z+1)+us3(z+1);
    %u2
         %Neacleated class
            un2(z+1) = sumun22+dt*(2*G(z)*sumun11);
         %Seed class
         %Virtual process
            us2(z+1) = sumus22+dt*(2*G(z)*sumus11);
         %the total crystal surface area
            uu2(z+1)=un2(z+1)+us2(z+1);

    %u1
         %Neacleated class
            un1(z+1) = sumun11+dt*(1*G(z)*sumun00);
         %Seed class
         %Virtual process
            us1(z+1) = sumus11+dt*(1*G(z)*sumus00);
         %the total crystal volume
            uu1(z+1)=un1(z+1)+us1(z+1);

    %u0
         %Neacleated class
            un0(z+1) = sumun00+dt*(0*G(z)*sumun00);
         %Seed class
         %Virtual process
            us0(z+1) = sumus00+dt*(0*G(z)*sumus00);
         %the total crystal volume
            uu0(z+1)=un0(z+1)+us0(z+1); 

    %u4
         %Neacleated class
            un4(z+1) = sumun44+dt*(4*G(z)*sumun44);
         %Seed class
         %Virtual process
            us4(z+1) = sumus44+dt*(4*G(z)*sumus44);
         %the total crystal volume
            uu4(z+1)=un4(z+1)+us4(z+1);

    %u5
         %Neacleated class
            un5(z+1) = sumun55+dt*(5*G(z)*sumun55);
         %Seed class
         %Virtual process
            us5(z+1) = sumus55+dt*(5*G(z)*sumus55);
         %the total crystal volume
            uu5(z+1)=un5(z+1)+us5(z+1);

end %fori=1:nt SS6=real(un3(p)); SS5=real(us3(p));

Lw=real(uu3/uu2); II1=[Tobj>=303]; II5= [Tobj<=323]; II2=[Cs<=Cobj]; II3= [us3(p)>=8.3301*10^9]; II4 = [Tobj(z+1)<=abs((2*dt)+Tobj(z))]

f= SS6/SS5;

   function [c1,c2,c3,c4,c5,ceq] = NONLCON(MV)

global Lw II1 II2 II3 II4 II5 Lw1 =Lw; c1 = double(II1) c2 = double(II2); c3 = double(II3); c3 = double(II4); c5 = double(II5) ceq = [];