Untitled

%Desarrollado por Juan Andrés Betancourt – Universidad Nacional de Colombia
%Ingeniería Química – Ingeniería de Procesos
%Algoritmo Genético de Optimización
%3 variables de optimización
%Cruce, mutación y selección elitista

%Tamaño de población
n = 80;
%Probabilidad operadores
cru = 0.9;
mut = 0.1;
%Generaciones
gen = 20;

%Rangos de variables
xinf = 0.29;
xsup = 0.33;
yinf = 0.05;
ysup = 0.15;
zinf = 25;
zsup = 35;

%Abrir Aspen Plus V10
Aspen = actxserver('Apwn.Document.36.0'); 
[stat,mess]=fileattrib; 
Simulation_Name = 'OPT';
Aspen.invoke('InitFromArchive2',[mess.Name '\' Simulation_Name '.bkp']);
Aspen.Visible = 0; 
Aspen.SuppressDialogs = 1; 

%Generar una población aleatoria de tamaño n
N = cell(n,1);
for k=1:length(N)
%Generar una matriz aleatoria 3x3
%Donde cada entrada es un entero aleatorio del 0 al 9
    I = zeros(3,3);
    for i=1:3
        for j=1:3
            I(i,j)=randi([0,9]);
        end
    end
    N{k}=I;
end
%Datos sobre unidades
%F3,5,8,12,24 en kg/s
%r3,4,5,12,14 en kg/m3
%QH1,2,3,4,5,6,7,REB,COND,WP1 en kW
%T6,8,7,11,23,24,16,17,29,30,31,32,18,33 en °C
%F28,22 en lb/h
%r28,22 en lb/ft3
%F19 en lb/min
%r19 en lb/gal
%BR en kmol/h
%F,F30CO2,NH3,F32CO2,NH3 en ton/h
%WW en fracción másica

%Inicio de algoritmo genético
d = 1
best=zeros(gen,1);
bestind=cell(gen,1);
while d<=gen
    %Realizar Cruce
    %Crear la población de hijos
    Nh = cell(n,1);
    for k=1:cru*n/2
    %Escoger dos individuos al azar:
        tag_I1 = randi([1,n]);
        tag_I2 = randi([1,n]);
        while tag_I1 == tag_I2
            tag_I2 = randi([1,n]);
        end
        I1 = N{tag_I1};
        I2 = N{tag_I2};
    %Escoger variable a cruzar
        col = randi([1,3]);
    %Escoger punto de cruce
        fil = randi([1,2]);
    %Cruzar I1 con I2
        I12 = I1;
        I12(fil+1:size(I12),col)=I2(fil+1:size(I12),col);
        I21 = I2;
        I21(fil+1:size(I21),col)=I1(fil+1:size(I21),col);
    %Agregar a población de hijos
        Nh{k}=I12;
        Nh{k+cru*n/2}=I21;
    end    
    %Realizar Mutación
    for k=1:mut*n
    %Escoger individuo al azar
        tag_I1 = randi([1,n]);
        I1 = N{tag_I1};
    %Escoger punto de mutación
        col = randi([1,3]);
        I11 = I1;
        if I11(size(I11,1),col)~=9
            I11(size(I11,1),col)=I11(size(I11,1),col)+1;
        else
            I11(size(I11,1),col)=0;
        end
    %Agregar a población de hijos
        Nh{k+cru*n}=I11;
    end
    %Evaluar la aptitud de padres e hijos
    M_APT = zeros(2*n,1);
    for k=1:n
    %Aptitud Padres
    x_=(N{k}(1,1)*100+N{k}(2,1)*10+N{k}(3,1))/1000;
    y_=(N{k}(1,2)*100+N{k}(2,2)*10+N{k}(3,2))/1000;
    z_=(N{k}(1,3)*100+N{k}(2,3)*10+N{k}(3,3))/1000;
    x = xinf+x_*(xsup-xinf);
    y = yinf+y_*(ysup-yinf);
    z = zinf+z_*(zsup-zinf);
    %Aspen.Reinit
    Aspen.Application.Tree.FindNode('\Data\Blocks\R-101\Input\CONV\2').value=x;
    Aspen.Application.Tree.FindNode('\Data\Blocks\S1\Input\FRAC\21').value=y;
    Aspen.Application.Tree.FindNode('\Data\Blocks\T-101\Input\T1').value=z;
    Run2(Aspen.Engine);
    F3=Aspen.Application.Tree.FindNode('\Data\Streams\3\Output\MASSFLMX\MIXED').value/3600;
    F5=Aspen.Application.Tree.FindNode('\Data\Streams\5\Output\MASSFLMX\MIXED').value/3600;
    F8=Aspen.Application.Tree.FindNode('\Data\Streams\8\Output\MASSFLMX\MIXED').value/3600;
    F12=Aspen.Application.Tree.FindNode('\Data\Streams\12\Output\MASSFLMX\MIXED').value/3600;
    F24=Aspen.Application.Tree.FindNode('\Data\Streams\24\Output\MASSFLMX\MIXED').value/3600;
    F28=Aspen.Application.Tree.FindNode('\Data\Streams\28\Output\MASSFLMX\MIXED').value*2.2046;
    F22=Aspen.Application.Tree.FindNode('\Data\Streams\22\Output\MASSFLMX\MIXED').value*2.2046;
    F19=Aspen.Application.Tree.FindNode('\Data\Streams\19\Output\MASSFLMX\MIXED').value*2.2046/60;
    r3=Aspen.Application.Tree.FindNode('\Data\Streams\3\Output\RHOMX_MASS\MIXED').value;
    r5=Aspen.Application.Tree.FindNode('\Data\Streams\5\Output\RHOMX_MASS\MIXED').value;
    r8=Aspen.Application.Tree.FindNode('\Data\Streams\8\Output\RHOMX_MASS\MIXED').value;
    r12=Aspen.Application.Tree.FindNode('\Data\Streams\12\Output\RHOMX_MASS\MIXED').value;
    r24=Aspen.Application.Tree.FindNode('\Data\Streams\24\Output\RHOMX_MASS\MIXED').value;
    r28=Aspen.Application.Tree.FindNode('\Data\Streams\28\Output\RHOMX_MASS\MIXED').value*2.2046/35.3147;
    r22=Aspen.Application.Tree.FindNode('\Data\Streams\22\Output\RHOMX_MASS\MIXED').value*2.2046/35.3147;
    r19=Aspen.Application.Tree.FindNode('\Data\Streams\19\Output\RHOMX_MASS\MIXED').value*2.2046/264.1720;
    QH1=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-101\Output\QCALC').value*4.184*10^6/3600;
    QH2=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-102\Output\QCALC').value*4.184*10^6/3600;
    QH3=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-103\Output\QCALC').value*4.184*10^6/3600;
    QH4=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-104\Output\QCALC').value*4.184*10^6/3600;
    QH5=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-105\Output\QCALC').value*4.184*10^6/3600;
    QH6=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-106\Output\QCALC').value*4.184*10^6/3600;
    QH7=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-107\Output\QCALC').value*4.184*10^6/3600;
    QREB=Aspen.Application.Tree.FindNode('\Data\Blocks\T-101\Output\REB_DUTY').value*4.184*10^6/3600;
    QCOND=Aspen.Application.Tree.FindNode('\Data\Blocks\T-101\Output\COND_DUTY').value*4.184*10^6/3600;
    WP1=Aspen.Application.Tree.FindNode('\Data\Blocks\T-101\Output\COND_DUTY').value;
    BR=Aspen.Application.Tree.FindNode('\Data\Blocks\T-101\Output\BOT_VFLOW').value;
    F=Aspen.Application.Tree.FindNode('\Data\Streams\27\Output\MASSFLOW\MIXED\UREA').value/1000;
    T6 = Aspen.Application.Tree.FindNode('\Data\Streams\6\Output\TEMP_OUT\MIXED').value;
    T8 = Aspen.Application.Tree.FindNode('\Data\Streams\8\Output\TEMP_OUT\MIXED').value;
    T7 = Aspen.Application.Tree.FindNode('\Data\Streams\7\Output\TEMP_OUT\MIXED').value;
    T11 = Aspen.Application.Tree.FindNode('\Data\Streams\11\Output\TEMP_OUT\MIXED').value;
    T23 = Aspen.Application.Tree.FindNode('\Data\Streams\23\Output\TEMP_OUT\MIXED').value;
    T24 = Aspen.Application.Tree.FindNode('\Data\Streams\24\Output\TEMP_OUT\MIXED').value;
    T16 = Aspen.Application.Tree.FindNode('\Data\Streams\16\Output\TEMP_OUT\MIXED').value;
    T17 = Aspen.Application.Tree.FindNode('\Data\Streams\17\Output\TEMP_OUT\MIXED').value;
    T29 = Aspen.Application.Tree.FindNode('\Data\Streams\29\Output\TEMP_OUT\MIXED').value;
    T30 = Aspen.Application.Tree.FindNode('\Data\Streams\30\Output\TEMP_OUT\MIXED').value;
    T31 = Aspen.Application.Tree.FindNode('\Data\Streams\31\Output\TEMP_OUT\MIXED').value;
    T32 = Aspen.Application.Tree.FindNode('\Data\Streams\32\Output\TEMP_OUT\MIXED').value;
    T18 = Aspen.Application.Tree.FindNode('\Data\Streams\18\Output\TEMP_OUT\MIXED').value;
    T33 = Aspen.Application.Tree.FindNode('\Data\Streams\33\Output\TEMP_OUT\MIXED').value;
    WW = Aspen.Application.Tree.FindNode('\Data\Streams\27\Output\MASSFRAC\MIXED\H2O').value;
    F30CO2 = Aspen.Application.Tree.FindNode('\Data\Streams\30\Output\MASSFLOW\MIXED\CO2').value/1000;
    F30NH3 = Aspen.Application.Tree.FindNode('\Data\Streams\30\Output\MASSFLOW\MIXED\NH3').value/1000;
    F32CO2 = Aspen.Application.Tree.FindNode('\Data\Streams\32\Output\MASSFLOW\MIXED\CO2').value/1000;
    F32NH3 = Aspen.Application.Tree.FindNode('\Data\Streams\32\Output\MASSFLOW\MIXED\NH3').value/1000;
    %Costos de Operación (USD/año)
    COHE1 = -6.69*QH1;
    COHE2 = 78.84*QH2;
    COHE3 = -6.69*QH3;
    COHE4 = 78.84*QH4;
    COHE5 = 78.84*QH5;
    COHE6 = 78.84*QH6;
    COHE7 = -6.69*QH7;
    COREB = 78.84*QREB;
    COCOND = -6.69*QCOND;
    COC1 = 0.034056*441.23*F28/r28;
    COC2 = 0.013651*441.23*F22/r22;
    COP1 = 678.9*WP1;
    CO = COHE1+COHE2+COHE3+COHE4+COHE5+COHE6+...
        COHE7+COREB+COCOND+COC1+COC2+COP1;
    %Costos de Capital (USD/año)
    CCR1 = 7474.576*(0.82*8*(4*((60*(-106.1921...
        *x^2+167.7313*x+9.2439))*F3/r3)/pi/8)...
        ^(1/3))^0.81*(1.093*(4*((60*(-106.1921*...
        x^2+167.7313*x+9.2439))*F3/r3)/pi/8)^...
        (1/3))^1.05;
    CCV1 = 7474.576*(0.82*2.5*(4*(300*F5/r5)/pi/2.5)...
        ^(1/3))^0.81*(1.093*(4*(300*F5/r5)/pi/2.5)^...
        (1/3))^1.05;
    CCV2 = 7474.576*(0.82*2.5*(4*(300*F8/r8)/pi/2.5)...
        ^(1/3))^0.81*(1.093*(4*(300*F8/r8)/pi/2.5)^...
        (1/3))^1.05;
    CCV3 = 7474.576*(0.82*2.5*(4*(300*F12/r12)/pi/2.5)...
        ^(1/3))^0.81*(1.093*(4*(300*F12/r12)/pi/2.5)^...
        (1/3))^1.05;
    CCD1 = 7474.576*(0.82*2.5*(4*(300*(F24/r24+400.15))/pi/2.5)...
        ^(1/3))^0.81*(1.093*(4*(300*(F24/r24+400.15))/pi/2.5)^...
        (1/3))^1.05;
    CCHE1 = 1387.1*(QH1/(2.555*(T8-T6))*10.7639)^0.65/5;
    CCHE2 = 1387.1*(QH2/(2.555*(T11-T7))*10.7639)^0.65/5;
    CCHE3 = 1007.1*(QH3/(2.555*(T24-T23))*10.7639)^0.65/5;
    CCHE4 = 1076.7*(QH4/(2.555*(T17-T16))*10.7639)^0.65/5;
    CCHE5 = 1663.1*(QH5/(2.555*(T30-T29))*10.7639)^0.65/5;
    CCHE6 = 1663.1*(QH6/(2.555*(T32-T31))*10.7639)^0.65/5;
    CCHE7 = 1076.7*(QH7/(2.555*(T33-T18))*10.7639)^0.65/5;
    CCREB = 1076.7*(QREB/(2.555*(250-T18))*10.7639)^0.65/5;
    CCCOND = 1076.7*(QCOND/(2.555*(-15))*10.7639)^0.65/5;
    CCC1 = 971.02*(0.034056*F28/r28)^0.82;
    CCC2 = 971.02*(0.013651*F22/r22)^0.82;
    CCP1 = (F19/r19)^0.64*549.705;
    CCT1 = 69313*(1.64*(BR*0.0001136+1.9597327))^1.45;
    CC = CCHE1+CCHE2+CCHE3+CCHE4+CCHE5+CCHE6+...
        CCHE7+CCREB+CCCOND+CCC1+CCC2+CCP1+CCT1+...
        CCR1+CCV1+CCV2+CCV3+CCD1;
    %Costos de MP
    CMPCO2 = (757.024-F30CO2-F32CO2)*24*365*176.6/0.975;
    CMPNH3 = (878.841-F30NH3-F32NH3)*24*365*262.5/0.995;
    CMP = CMPCO2+CMPNH3;
    %Ingresos por venta de urea (USD/año) 
    Iv = F*477.5*24*365; 
    %F. O. (MUSD/año)
    FO = (Iv-CC-CO-CMP)/10^6;
    rest = 0.005-WW;
    if rest >= 0
        APT = FO;
    else 
        APT = FO*0.8;
    end
    M_APT(k,1)=APT;
    %Aptitud Hijos
    x_=(Nh{k}(1,1)*100+Nh{k}(2,1)*10+Nh{k}(3,1))/1000;
    y_=(Nh{k}(1,2)*100+Nh{k}(2,2)*10+Nh{k}(3,2))/1000;
    z_=(Nh{k}(1,3)*100+Nh{k}(2,3)*10+Nh{k}(3,3))/1000;
    x = xinf+x_*(xsup-xinf);
    y = yinf+y_*(ysup-yinf);
    z = zinf+z_*(zsup-zinf);
    %Aspen.Reinit
    Aspen.Application.Tree.FindNode('\Data\Blocks\R-101\Input\CONV\2').value=x;
    Aspen.Application.Tree.FindNode('\Data\Blocks\S1\Input\FRAC\21').value=y;
    Aspen.Application.Tree.FindNode('\Data\Blocks\T-101\Input\T1').value=z;
    Run2(Aspen.Engine);
    F3=Aspen.Application.Tree.FindNode('\Data\Streams\3\Output\MASSFLMX\MIXED').value/3600;
    F5=Aspen.Application.Tree.FindNode('\Data\Streams\5\Output\MASSFLMX\MIXED').value/3600;
    F8=Aspen.Application.Tree.FindNode('\Data\Streams\8\Output\MASSFLMX\MIXED').value/3600;
    F12=Aspen.Application.Tree.FindNode('\Data\Streams\12\Output\MASSFLMX\MIXED').value/3600;
    F24=Aspen.Application.Tree.FindNode('\Data\Streams\24\Output\MASSFLMX\MIXED').value/3600;
    F28=Aspen.Application.Tree.FindNode('\Data\Streams\28\Output\MASSFLMX\MIXED').value*2.2046;
    F22=Aspen.Application.Tree.FindNode('\Data\Streams\22\Output\MASSFLMX\MIXED').value*2.2046;
    F19=Aspen.Application.Tree.FindNode('\Data\Streams\19\Output\MASSFLMX\MIXED').value*2.2046/60;
    r3=Aspen.Application.Tree.FindNode('\Data\Streams\3\Output\RHOMX_MASS\MIXED').value;
    r5=Aspen.Application.Tree.FindNode('\Data\Streams\5\Output\RHOMX_MASS\MIXED').value;
    r8=Aspen.Application.Tree.FindNode('\Data\Streams\8\Output\RHOMX_MASS\MIXED').value;
    r12=Aspen.Application.Tree.FindNode('\Data\Streams\12\Output\RHOMX_MASS\MIXED').value;
    r24=Aspen.Application.Tree.FindNode('\Data\Streams\24\Output\RHOMX_MASS\MIXED').value;
    r28=Aspen.Application.Tree.FindNode('\Data\Streams\28\Output\RHOMX_MASS\MIXED').value*2.2046/35.3147;
    r22=Aspen.Application.Tree.FindNode('\Data\Streams\22\Output\RHOMX_MASS\MIXED').value*2.2046/35.3147;
    r19=Aspen.Application.Tree.FindNode('\Data\Streams\19\Output\RHOMX_MASS\MIXED').value*2.2046/264.1720;
    QH1=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-101\Output\QCALC').value*4.184*10^6/3600;
    QH2=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-102\Output\QCALC').value*4.184*10^6/3600;
    QH3=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-103\Output\QCALC').value*4.184*10^6/3600;
    QH4=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-104\Output\QCALC').value*4.184*10^6/3600;
    QH5=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-105\Output\QCALC').value*4.184*10^6/3600;
    QH6=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-106\Output\QCALC').value*4.184*10^6/3600;
    QH7=Aspen.Application.Tree.FindNode('\Data\Blocks\HE-107\Output\QCALC').value*4.184*10^6/3600;
    QREB=Aspen.Application.Tree.FindNode('\Data\Blocks\T-101\Output\REB_DUTY').value*4.184*10^6/3600;
    QCOND=Aspen.Application.Tree.FindNode('\Data\Blocks\T-101\Output\COND_DUTY').value*4.184*10^6/3600;
    WP1=Aspen.Application.Tree.FindNode('\Data\Blocks\T-101\Output\COND_DUTY').value;
    BR=Aspen.Application.Tree.FindNode('\Data\Blocks\T-101\Output\BOT_VFLOW').value;
    F=Aspen.Application.Tree.FindNode('\Data\Streams\27\Output\MASSFLOW\MIXED\UREA').value/1000;
    T6 = Aspen.Application.Tree.FindNode('\Data\Streams\6\Output\TEMP_OUT\MIXED').value;
    T8 = Aspen.Application.Tree.FindNode('\Data\Streams\8\Output\TEMP_OUT\MIXED').value;
    T7 = Aspen.Application.Tree.FindNode('\Data\Streams\7\Output\TEMP_OUT\MIXED').value;
    T11 = Aspen.Application.Tree.FindNode('\Data\Streams\11\Output\TEMP_OUT\MIXED').value;
    T23 = Aspen.Application.Tree.FindNode('\Data\Streams\23\Output\TEMP_OUT\MIXED').value;
    T24 = Aspen.Application.Tree.FindNode('\Data\Streams\24\Output\TEMP_OUT\MIXED').value;
    T16 = Aspen.Application.Tree.FindNode('\Data\Streams\16\Output\TEMP_OUT\MIXED').value;
    T17 = Aspen.Application.Tree.FindNode('\Data\Streams\17\Output\TEMP_OUT\MIXED').value;
    T29 = Aspen.Application.Tree.FindNode('\Data\Streams\29\Output\TEMP_OUT\MIXED').value;
    T30 = Aspen.Application.Tree.FindNode('\Data\Streams\30\Output\TEMP_OUT\MIXED').value;
    T31 = Aspen.Application.Tree.FindNode('\Data\Streams\31\Output\TEMP_OUT\MIXED').value;
    T32 = Aspen.Application.Tree.FindNode('\Data\Streams\32\Output\TEMP_OUT\MIXED').value;
    T18 = Aspen.Application.Tree.FindNode('\Data\Streams\18\Output\TEMP_OUT\MIXED').value;
    T33 = Aspen.Application.Tree.FindNode('\Data\Streams\33\Output\TEMP_OUT\MIXED').value;
    WW = Aspen.Application.Tree.FindNode('\Data\Streams\27\Output\MASSFRAC\MIXED\H2O').value;
    F30CO2 = Aspen.Application.Tree.FindNode('\Data\Streams\30\Output\MASSFLOW\MIXED\CO2').value/1000;
    F30NH3 = Aspen.Application.Tree.FindNode('\Data\Streams\30\Output\MASSFLOW\MIXED\NH3').value/1000;
    F32CO2 = Aspen.Application.Tree.FindNode('\Data\Streams\32\Output\MASSFLOW\MIXED\CO2').value/1000;
    F32NH3 = Aspen.Application.Tree.FindNode('\Data\Streams\32\Output\MASSFLOW\MIXED\NH3').value/1000;
    %Costos de Operación (USD/año)
    COHE1 = -6.69*QH1;
    COHE2 = 78.84*QH2;
    COHE3 = -6.69*QH3;
    COHE4 = 78.84*QH4;
    COHE5 = 78.84*QH5;
    COHE6 = 78.84*QH6;
    COHE7 = -6.69*QH7;
    COREB = 78.84*QREB;
    COCOND = -6.69*QCOND;
    COC1 = 0.034056*441.23*F28/r28;
    COC2 = 0.013651*441.23*F22/r22;
    COP1 = 678.9*WP1;
    CO = COHE1+COHE2+COHE3+COHE4+COHE5+COHE6+...
        COHE7+COREB+COCOND+COC1+COC2+COP1;
    %Costos de Capital (USD/año)
    CCR1 = 7474.576*(0.82*8*(4*((60*(-106.1921...
        *x^2+167.7313*x+9.2439))*F3/r3)/pi/8)...
        ^(1/3))^0.81*(1.093*(4*((60*(-106.1921*...
        x^2+167.7313*x+9.2439))*F3/r3)/pi/8)^...
        (1/3))^1.05;
    CCV1 = 7474.576*(0.82*2.5*(4*(300*F5/r5)/pi/2.5)...
        ^(1/3))^0.81*(1.093*(4*(300*F5/r5)/pi/2.5)^...
        (1/3))^1.05;
    CCV2 = 7474.576*(0.82*2.5*(4*(300*F8/r8)/pi/2.5)...
        ^(1/3))^0.81*(1.093*(4*(300*F8/r8)/pi/2.5)^...
        (1/3))^1.05;
    CCV3 = 7474.576*(0.82*2.5*(4*(300*F12/r12)/pi/2.5)...
        ^(1/3))^0.81*(1.093*(4*(300*F12/r12)/pi/2.5)^...
        (1/3))^1.05;
    CCD1 = 7474.576*(0.82*2.5*(4*(300*(F24/r24+400.15))/pi/2.5)...
        ^(1/3))^0.81*(1.093*(4*(300*(F24/r24+400.15))/pi/2.5)^...
        (1/3))^1.05;
    CCHE1 = 1387.1*(QH1/(2.555*(T8-T6))*10.7639)^0.65/5;
    CCHE2 = 1387.1*(QH2/(2.555*(T11-T7))*10.7639)^0.65/5;
    CCHE3 = 1007.1*(QH3/(2.555*(T24-T23))*10.7639)^0.65/5;
    CCHE4 = 1076.7*(QH4/(2.555*(T17-T16))*10.7639)^0.65/5;
    CCHE5 = 1663.1*(QH5/(2.555*(T30-T29))*10.7639)^0.65/5;
    CCHE6 = 1663.1*(QH6/(2.555*(T32-T31))*10.7639)^0.65/5;
    CCHE7 = 1076.7*(QH7/(2.555*(T33-T18))*10.7639)^0.65/5;
    CCREB = 1076.7*(QREB/(2.555*(250-T18))*10.7639)^0.65/5;
    CCCOND = 1076.7*(QCOND/(2.555*(-15))*10.7639)^0.65/5;
    CCC1 = 971.02*(0.034056*F28/r28)^0.82;
    CCC2 = 971.02*(0.013651*F22/r22)^0.82;
    CCP1 = (F19/r19)^0.64*549.705;
    CCT1 = 69313*(1.64*(BR*0.0001136+1.9597327))^1.45;
    CC = CCHE1+CCHE2+CCHE3+CCHE4+CCHE5+CCHE6+...
        CCHE7+CCREB+CCCOND+CCC1+CCC2+CCP1+CCT1+...
        CCR1+CCV1+CCV2+CCV3+CCD1;
    %Costos de MP
    CMPCO2 = (757.024-F30CO2-F32CO2)*24*365*176.6/0.975;
    CMPNH3 = (878.841-F30NH3-F32NH3)*24*365*262.5/0.995;
    CMP = CMPCO2+CMPNH3;
    %Ingresos por venta de urea (USD/año) 
    Iv = F*477.5*24*365; 
    %F. O. (MUSD/año)
    FO = (Iv-CC-CO-CMP)/10^6;
    rest = 0.005-WW;
    if rest >= 0
        APT = FO;
    else 
        APT = FO*0.8;
    end
    M_APT(k+n)=APT;
    end
    %Selección elitista
    %Ordenar aptitudes y conservar índice
    [Mcut,Index] = sort(M_APT);
    %Generar la nueva población de padres
    Nnew = cell(n,1);
    %Asignar los nuevos padres según aptitud
    l=1;
    for k=n+1:2*n
    Ind = Index(k);
    if Ind > n
        Nnew{l}=Nh{Ind-n};
    else
        Nnew{l}=N{Ind};
    end
    l=l+1;
    end
    best(d) = Mcut(2*n);
    Ibest = Index(2*n);
    if Ibest > n
        bestind{d}=Nh{Ibest-n};
    else
        bestind{d}=N{Ibest};
    end
    bestind{d}
    best(d)
    N = Nnew;
    d = d+1
end

%Traducir los datos de la corrida
individuos= zeros(gen,5); 
for k=1:gen
    x_=(bestind{k}(1,1)*100+bestind{k}(2,1)*10+bestind{k}(3,1))/1000;
    y_=(bestind{k}(1,2)*100+bestind{k}(2,2)*10+bestind{k}(3,2))/1000;
    z_=(bestind{k}(1,3)*100+bestind{k}(2,3)*10+bestind{k}(3,3))/1000;
    x = xinf+x_*(xsup-xinf);
    y = yinf+y_*(ysup-yinf);
    z = zinf+z_*(zsup-zinf);
    individuos(k,1)=k;
    individuos(k,2)=best(k);
    individuos(k,3)=x;
    individuos(k,4)=y;
    individuos(k,5)=z;
end