hookes pattern search optimisaton

Question

sai.aleenus 0 Newbie Poster

11 Years Ago

can someone help me with the development of a hook and jeeves pattern search optimisation in c programming with graphical visualisation in step wise.
a function bhp is to be minimised by effective selection of 2 parameters x1 and x2 while bhp has been passed x1 ans x2 as arguments of the func.

c

2 Contributors
32 Replies
693 Views
1 Week Discussion Span
Latest Post 11 Years Ago Latest Post by Schol-R-LEA

All 32 Replies

Schol-R-LEA 1,446 Commie Mutant Traitor

11 Years Ago

I would recommend moving the initialization data into a file to be loaded at runtime; this would allow you to have different starting values for the plant's state variables, without having to alter the compiled program. You can have different files containing different starting values to get the variance you need.

While I dislike having so many global variables, I don't see any alternative, as there is no clear way to pass so many values to the function without either an unreasonably large number of parameters, or lumping them all together in a massive struct. As it is, the best we can hope for is to rein in the chaos by putting related values together into arrays and structs, and even then, I don't know enough about the data sets to give meaningful names and develop useable groupings to all the data.

I have managed to pull together some other functions, and at least one useable structure, but I don't know how helpful these will be to you. Note that I have been using double variables rather than floats, to get higher precision; I would recommend you do the same throughout the program.

struct
{
    double hum1, hum2, po, c;
} humidity;


double sp_humidity(double d, double w, double svp)
{
    double sphua, sphub, sphud;

    sphua = humidity.hum1 / humidity.hum2;
    sphub = humidity.po * humidity.c * (d - w);
    sphud = svp - sphub;
    return (sphua * sphud) / (humidity.po - sphud); // sp humidity is calculated
}


double enthalpy_metric(double d, double humidity)
{
    return (spec2 * d) + (humidity * (spelo + (spec3 * d)));
}

double enthalpy_btu(double enthalpy)
{
    return ((enthalpy * 1.8) / 4.1858) + 7.686;
}

On a related note, these all seem to be called on a series of values named things like spe_in, spe1, spe2, etc. such that they form groups of five values. This can be simplified thusly:

    double svp[5];
    double hin_f[5];
    double twbi[5], svp[5], t[5];
    double sphu[5],spe[5];


    twbi[0] = 37.0f;
    twbi[1] = twbi[0] + (0.1*r);
    twbi[2] = twbi[0] + (0.4*r);
    twbi[3] = twbi[0] + (0.6*r);
    twbi[4] = twbi[0] + (0.9*r);

    /*   then later in the code ... */

    for (i = 0; i < 5; i++)
    {
        svp[i] = pow(10, saturation_vapor_pressure(twbi[i]));                //sat vap pressure calculated
        printf("\n sat vap pressure is %f Pa",svp[i]);

        sphu[i] = sp_humidity(tdbi, twbi[i], svp[i]);
        printf("\n sp humidity is %f kg/kg of dry air",sphu[i]);

        spe[i] = enthalpy_metric(tdbi, sphu[i]);//sp enthalpy is calculated
        printf("\n enthalpy is(kj/kg) %f",spe[i]);

        hin_f[i] = enthalpy_btu(spe[i]);
        printf("\n enthalpy is(btu/lb) %f", hin_f[i]);
    }

This moves the entire group out of the larger loop, and forms a logical grouping of data points.

If you don't mind me asking, what is your major? For some reason, I get the impression you are a chemical or industrial engineering student rather than a CS major per se.

Edited 11 Years Ago by Schol-R-LEA

Schol-R-LEA 1,446 Commie Mutant Traitor

11 Years Ago

OK, now, I have done some major surgery on the bhpeqns() function - purely refactoring, mind you, it should give more or less the same results as before - which should make fixing and maintaining the code much, much easier. In the process, I've managed to cut the number of local variables in the function by about 2/3rds. If I actually understood the code better, I could probably do more still. :)

The two major realizations which made this possible were the observations that several of the variables were for the same values, but in metric and US units, of which many were for either the intake valve or the outflow valve. By defining a struct type, I was able to combine the first group into single variables; by changing the second ones to arrays, with defined indices for the INLET and OUTLET; I was able to combine the latter. In many cases, both applied, so I was able to replace four separate variables with a single composite variable.

Note that I systematically used double rather than float, as you get significant improvements in the calculations' precision this way.

You will also be pleasantly surprised by the way I have simplified the calculations for lg, KAV1 and KAV2, I think.

I think that you would do well to adopt this version of the function, as it should be much easier to understand once you've grasped the way I've organized the variables, and much easier to fix any problems that may arise in it. I've checked that it gives the same results for at least those datasets I've tested so far.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdbool.h>
#include "hooke-jeeves.h"


#define VARS (250)
#define RHO (0.85)
#define EPSILON (0.0001)
#define ITERMAX (1024)
#define NVARS (2)


#define INLET 0
#define OUTLET 5

struct DUAL_MEASUREMENT
{
    double metric, us;
};

bool initialize_plant_state(FILE* data_src);

double bhpeqns(double fh, double fs);

double saturation_vapor_pressure(double twbi, bool at_exit);
double sp_humidity(double d, double w, double svp);
double sp_volume(double tdb, double twb, double svp);
double enthalpy_metric(double d, double humid);
double enthalpy_btu(double enthalpy);
double kgs_to_lbsmin(double flow_rate);
double cums_to_cfts(double cums);
double compute_kav1(double lg, double r, struct DUAL_MEASUREMENT spe[], double r_offset[]);
double compute_kav2(double lg, double fvel, double fh, double fs);


// global 'constants' to be set in initialize_plant_state() -
// they vary with the initial data set, but are not altered during the
// program run

double thw, tcw, twb_init, tdbi, fand, hubd, fann;
double aih, cell, celw, gpm;
int ais, celn;
double to, t9;
double vapor[8];
struct
{
    double hum1, hum2, po, c;
} humidity;

double spec[3], spvolro;
int wlmin, wlmax, vmin, vmax;


// global variables - these are set in the functions
double app;
struct DUAL_MEASUREMENT edisa, wl;
double aia, cela, gpmpc;


int main()
{
    char filename[FILENAME_MAX];
    FILE* data_src;
    bool success;

    printf("Enter the name of the data file: ");
    scanf("%4095s", filename);
    data_src = fopen(filename, "r");

    if (data_src == NULL)
    {
        printf("Unable to open file '%s', exiting.\n", filename);
        exit(1);
    }

    success = initialize_plant_state(data_src);
    fclose(data_src);

    if (!success)
    {
        printf("Could not read all the require data, exiting.\n");
        exit(1);
    }

    bhpeqns(0.6096, 12.0);

    return 0;
}



bool initialize_plant_state(FILE* data_src)
{
    if (7 > fscanf(data_src, "%lf %lf %lf %lf %lf %lf %lf", &thw, &tcw, &twb_init, &tdbi, &fand, &hubd, &fann))
    {
        printf("line 0: could not read one or more of the following: %lf %lf %lf %lf %lf %lf %lf\n", thw, tcw, twb_init, tdbi, fand, hubd, fann);
        return false;
    }

    if (4 > fscanf(data_src, "%lf %lf %lf %lf\n", &aih, &cell, &celw, &gpm))
    {
        printf("line 1: could not read one or more of the following: %lf %lf %lf %lf\n", aih, cell, celw, gpm);
        return false;
    }
    if (2 > fscanf(data_src, "%d %d\n", &ais, &celn))
    {
        printf("line 2: could not read one or more of the following: %d %d\n", ais, celn);
        return false;
    }
    if (2 > fscanf(data_src, "%lf %lf\n", &to, &t9))
    {
        printf("line 3: could not read one or more of the following: %lf %lf\n", to, t9);
        return false;
    }
    if (7 > fscanf(data_src, "%lf %lf %lf %lf %lf %lf %lf %lf\n", &vapor[0], &vapor[1], &vapor[2], &vapor[3], &vapor[4], &vapor[5], &vapor[6], &vapor[7]))
    {
        printf("line 4: could not read one or more of the following: %lf %lf %lf %lf %lf %lf %lf %lf\n",
               vapor[0], vapor[1], vapor[2], vapor[3], vapor[4], vapor[5], vapor[6], vapor[7]);
        return false;
    }
    if (4 > fscanf(data_src, "%lf %lf %lf %lf\n", &humidity.hum1, &humidity.hum2, &humidity.po, &humidity.c))
    {
        printf("line 5: could not read one or more of the following: %lf %lf %lf %lf\n", humidity.hum1, humidity.hum2, humidity.po, humidity.c);
        return false;
    }
    if (4 > fscanf(data_src, "%lf %lf %lf %lf\n", &spec[1], &spec[2], &spec[0], &spvolro))
    {
        printf("line 6: could not read one or more of the following: %lf %lf %lf %lf\n", spec[1], spec[2], spec[0], spvolro);
        return false;
    }
    if (4 > fscanf(data_src, "%d %d %d %d\n", &wlmin, &wlmax, &vmin, &vmax))
    {
        printf("line 7: could not read one or more of the following: %d %d %d %d\n", wlmin, wlmax, vmin, vmax);
        return false;
    }

    return true;
}


double saturation_vapor_pressure(double twbi, bool at_exit)
{
    double too, saturation[8];

    too = twbi + to;
    saturation[0] = vapor[0] / too;
    saturation[1] = vapor[1] * log(too);
    saturation[2] = vapor[3] * (twbi - 0.01);
    saturation[3] = vapor[2] * pow(10, saturation[2]);
    saturation[4] = vapor[5] * (1 - (t9 / too));
    saturation[5] = vapor[4] * pow(10, saturation[4]);
    saturation[6] = vapor[6] + saturation[0] + saturation[1] + saturation[3] + saturation[5];

    if (at_exit)
    {
        saturation[7] = saturation[6] + vapor[7];
    }
    else
    {
        saturation[7] = saturation[6];
    }

    return pow(10, saturation[7]);                //sat vap pressure calculated
}



double sp_humidity(double d, double w, double svp)
{
    double hum[3];

    hum[0] = humidity.hum1 / humidity.hum2;
    hum[1] = humidity.po * humidity.c * (d - w);
    hum[2] = svp - hum[1];
    return fabs((hum[0]* hum[2]) / (humidity.po - hum[2])); // sp humidity is calculated
}


double sp_volume(double tdb, double twb, double svp)
{
    const double sphup[] = {18.01534, 28.9645,  101325,  0.000666};
    double vol[5],fspvol;

    vol[0] = spvolro * (tdb + to);
    vol[1] = sphup[2] * sphup[3] * (tdb - twb);
    vol[2] = sphup[2] - svp + vol[1];
    vol[3] = sphup[1] * vol[2] ;
    vol[4] = (vol[0] * 1000)/vol[3] ;
    fspvol = vol[4];
    fspvol = fabs(fspvol);
    return(fspvol);
}


double enthalpy_metric(double d, double humid)
{
    return (spec[1] * d) + (humid * (spec[0] + (spec[2] * d)));
}

double enthalpy_btu(double enthalpy)
{
    return ((enthalpy * 1.8) / 4.1858) + 7.686;
}

double kgs_to_lbsmin(double flow_rate)
{
    return flow_rate / 0.00755987283;
}

double cums_to_cfts(double cums)
{
    return cums * pow(3.2808, 2);
}

double compute_kav1(double lg, double r, struct DUAL_MEASUREMENT spe[], double r_offset[])
{
    int i;
    double speat, ra, ha, hb, h, h_inv, h_inv_sum;

    for (i = 0; i < 4; i++)
    {
        speat = spe[INLET].metric + (r_offset[i] * lg * r);
        ha = spe[i + 1].metric - speat;
        hb = ((ha * 1.8) / 4.1858);
        h = hb + 7.686 ;
        h_inv = 1 / h;
        h_inv_sum += h_inv;
    }

    ra = (r * 1.8) + 32;//in deg f

    return ra * 0.25 * h_inv_sum;   //in fps ---- demand
}


double compute_kav2(double lg, double fvel, double fh, double fs)
{
    double kav2_sub;

    kav2_sub = 0.06717927 * pow(lg, -0.546673);
    kav2_sub *= pow(fvel, 0.25158684);
    kav2_sub *= pow(fh, 0.6369974);
    kav2_sub *= pow (fs, -0.429);
    return kav2_sub;
}


double bhpeqns(double fh_init, double fs_init)
{
    int i;

    double r, r_offset[] = {0.1, 0.4, 0.6, 0.9};

    struct DUAL_MEASUREMENT fh, fs;
    double twb[OUTLET + 1];
    double svp[OUTLET + 1], sphu[OUTLET + 1];
    struct DUAL_MEASUREMENT spe[OUTLET + 1];
    double kav1, kav2, lg;
    struct DUAL_MEASUREMENT fvel, liquid_flow_rate, gas_flow_rate;
    struct DUAL_MEASUREMENT spvol[OUTLET + 1], row[OUTLET + 1], vel[OUTLET + 1];
    struct DUAL_MEASUREMENT spvol_fill, row_fill, vel_fill;
    double pdf[6];
    struct DUAL_MEASUREMENT cfm_fan, pd_v, tpd, bhp;
    double pd_il_inh20, pd_de_inh20, pd_fill_inh20;
    double tpd_s_pa, tpd_s_inh20;
    double spvolref, spvolre;

    fh.metric = fh_init;
    fs.metric = fs_init;

    r = thw - tcw;
    app = tcw - twb[INLET] ;
    cela = cell * celw;
    aia = ais * celw * aih;
    edisa.metric = 0.785 * (pow(fand, 2) - pow(hubd, 2));
    edisa.us = cums_to_cfts(edisa.metric);
    printf("\n edisa %lf", edisa.us);
    gpmpc = gpm / celn;
    wl.metric = gpmpc / cela;
    wl.us = (gpmpc * 4.402867) / cums_to_cfts(cela);  // fpm

    twb[INLET] = twb_init;
    twb[1] = twb[INLET] + (r_offset[0] * r);
    twb[2] = twb[INLET] + (r_offset[1] * r);
    twb[3] = twb[INLET] + (r_offset[2] * r);
    twb[4] = twb[INLET] + (r_offset[3] * r);


    for (i = 0; i < OUTLET; i++)
    {

        svp[i] = saturation_vapor_pressure(twb[i], false);
        sphu[i] = sp_humidity(tdbi, twb[i], svp[i]);
        spe[i].metric = enthalpy_metric(tdbi, sphu[i]);//sp enthalpy is calculated
        spe[i].us = enthalpy_btu(spe[i].metric);

        if (i == INLET)
        {
            printf("\n inlet sat vap pressure is %lf Pa", svp[INLET]);
            printf("\n inlet sp humidity is %lf kg/kg of dry air", sphu[INLET]);
            printf("\n inlet enthalpy is(kj/kg) %lf", spe[INLET].metric);
            printf("\n inlet enthalpy is(btu/lb) %lf", spe[INLET].us);
        }
    }

    for(lg = 1.0, kav1 = 0, kav2 = 1; ((kav1 <= kav2) && (lg <= 3) && (lg >= 1)); lg += 0.001)
    {
        kav1 = compute_kav1(lg, r, spe, r_offset);

        fvel.us = (wl.metric * 8.33) / (lg * 0.07 * 0.227125);  //in fpm conversion factor 0.227125 used
        fh.us = fh.metric / 0.3048;
        fs.us = fs.metric * 0.00328084;  //in ft
        kav2 = compute_kav2(lg, fvel.us, fh.us, fs.us);
    }


    printf("\n the l/g ratio required for the performance is %lf \n", lg);
    printf("\n the KAV/L_2  for the performance is %lf \n", kav2);
    printf("\n the KAV/L_1  for the performance is %lf \n", kav1);


    fvel.us = (wl.metric * 8.33) / (lg * 0.07 * 0.00508);

    liquid_flow_rate.metric = (gpmpc * 1000) / 3600;// kg/sec
    liquid_flow_rate.us = kgs_to_lbsmin(liquid_flow_rate.metric);
    gas_flow_rate.metric = liquid_flow_rate.metric / lg;
    gas_flow_rate.us = kgs_to_lbsmin(gas_flow_rate.metric);
    printf("\n the mass flow rate of water is %lf kg/s or %lf lb/min \n", liquid_flow_rate.metric, liquid_flow_rate.us);

    printf("\n the mass flow rate of air is %lf kg/s or %lf lb/min \n", gas_flow_rate.metric, gas_flow_rate.us);


    printf("\n AT COOLING TOWER INLET\n");

    printf("\n sat vap p inlet is %lf Pa", svp[0]);
    printf("\n sp humidity inlet is %lf kg/kg of dry air", sphu[0]);
    spvol[INLET].metric = sp_volume(tdbi, twb[INLET], svp[INLET]);
    spvol[INLET].us = spvol[INLET].metric / 0.0624279606; //cuft/lb
    printf("\n the inlet specific volume is %lf cum/kg or %lf cu ft/lb ",spvol[INLET].metric, spvol[INLET].us);//

    row[INLET].metric =(1+ sphu[0])/(spvol[INLET].metric);
    row[INLET].us = row[INLET].metric / 16.0184634;
    printf("\n the inlet density is %lf kg/cum or %lf lb/cu ft", row[INLET].metric, row[INLET].us);


    spe[OUTLET].metric = spe[INLET].metric + (lg * r);
    spe[OUTLET].us = enthalpy_btu(spe[OUTLET].metric);
    printf("\n the inlet enthalpy is %lf kj/kg or %lf btu/lb", spe[INLET].metric, spe[INLET].us);
    printf("\n the exit enthalpy is %lf kj/kg or %lf btu/lb", spe[OUTLET].metric, spe[OUTLET].us);

    twb[OUTLET] = spe[OUTLET].metric / 1.006;
    printf("\n twb of the outlet is %lf", twb[OUTLET]);

    svp[OUTLET] = saturation_vapor_pressure(twb[OUTLET], true);
    printf("\n sat vap p exit is %lf Pa   svpn    %lf", svp[OUTLET], svp[INLET]);

    sphu[OUTLET] = sp_humidity(twb[OUTLET], twb[OUTLET], svp[OUTLET]);
    printf("\n sp humidity exit is %lf kg/kg of saturated air", sphu[OUTLET]);


    printf("\n AT COOLING TOWER EXIT\n");

    spvolre = svp[OUTLET] - humidity.po;
    spvolref = abs(spvolre);
    spvol[OUTLET].metric = (spvolro * (twb[OUTLET] + to) * 1000) / (humidity.hum2 * spvolref);  //cum/kg
    spvol[OUTLET].us = spvol[OUTLET].metric / 0.0624279606;//cuft/lb
    printf("\n the exit specific volume is %lf cum/kg or %lf cu ft/lb ",spvol[OUTLET].metric, spvol[OUTLET].us);

    row[OUTLET].metric =(1+ sphu[OUTLET])/(spvol[OUTLET].metric);

    row[OUTLET].us = row[OUTLET].metric/16.0184634;
    printf("\n the exit density is %lf kg/cum or %lf lb/cu ft",row[OUTLET].metric,row[OUTLET].us);

    vel[INLET].metric = (gas_flow_rate.metric * spvol[INLET].metric) / aia ;//m/s
    vel[INLET].us = vel[INLET].metric/0.00508;//fpm
    vel[OUTLET].metric = (gas_flow_rate.metric * spvol[OUTLET].metric) / cela ;
    vel[OUTLET].us =  vel[OUTLET].metric / 0.00508;

    row_fill.metric = (row[INLET].metric + row[OUTLET].metric) / 2;
    row_fill.us = (row[INLET].us + row[OUTLET].us)/2 ;
    spvol_fill.metric = (spvol[INLET].metric + spvol[OUTLET].metric) / 2;
    spvol_fill.us = (spvol[INLET].us + spvol[OUTLET].us) / 2;   //

    vel_fill.metric = (gas_flow_rate.metric * spvol_fill.metric) / cela;
    vel_fill.us = vel_fill.metric / 0.00508;

    printf("\n the specific volume at fill is %lf cum/kg or %lf cu ft/lb ", spvol_fill.metric, spvol_fill.us);
    printf("\n the density at fill is %lf kg/cum or %lf lb/cu ft", row_fill.metric, row_fill.us);
    printf("\n the velocity of air at inlet is %lf m/s or %lf fpm ", vel[INLET].metric, vel[INLET].us);
    printf("\n the velocity of air at fill is %lf m/s or %lf fpm ", vel_fill.metric, vel_fill.us);
    printf("\n the velocity of air at exit is %lf m/s or %lf fpm ", vel[OUTLET].metric, vel[OUTLET].us);

    cfm_fan.metric = gas_flow_rate.metric * spvol[OUTLET].metric;//cum/s
    cfm_fan.us = gas_flow_rate.us * spvol[OUTLET].us;
    printf("\n the volume flow rate of air at fan CFAN is %lf cum/s or %lf cu ft/min \n", cfm_fan.metric, cfm_fan.us);



    pd_il_inh20 = (3 * ((2 * pow(10,-7) * vel[INLET].us * vel[INLET].us) + (3 * pow(10,-6) * vel[INLET].us)) * (row[INLET].us / 0.07));
    //       pd_il_pa = pd_il_inh20*249;//pascal
    printf("\n pd il %lf in ", pd_il_inh20);

    pd_de_inh20 = (5 * ((3 * pow(10,-7) * vel[OUTLET].us * vel[OUTLET].us) - (9 * pow(10,-6) * vel[OUTLET].us) + 0.003 ) * (row[OUTLET].us / 0.07));
    //       pd_de_pa = pd_de_inh20*249;//pascal
    printf("\n pd de %lf in ", pd_de_inh20);



    pdf[0] = (fh.us / fs.us);
    pdf[1] = pow(pdf[0], 0.963077);
    pdf[2] = ((7.44385 * pow(10,-6)) + ((-4.078 * pow(10,-6)) * pow(fh.us, 0.176686)));
    pdf[3] = (4.897 * (pow(10,-4)) * (pow(vel_fill.us, 2.3386)));
    pdf[4] = (0.0428 * wl.us);
    pdf[5] = (pow(vel_fill.us, 1.33558));

    pd_fill_inh20 = (pdf[1] * pdf[2] * (pdf[3] + (pdf[4] * pdf[5]))) * (row_fill.us / 0.07);
    printf("\n pd fill %lf in ", pd_fill_inh20);
    //      pd_fill_pa = pd_fill_inh20*249;

    tpd_s_inh20 = pd_il_inh20 + pd_de_inh20 + pd_fill_inh20;
    printf("\n hs %lf in ", tpd_s_inh20);
    tpd_s_pa = tpd_s_inh20 * 249;
    pd_v.us = (row[OUTLET].us / 1202312) * (cfm_fan.us / edisa.us) * (cfm_fan.us / edisa.us);
    printf("\n hv %lf in", pd_v.us);
    pd_v.metric = pd_v.us * 249 ;
    tpd.us = pd_v.us + tpd_s_inh20;
    tpd.metric = tpd_s_pa + pd_v.metric;
    printf("\n the overall pressure drop across the cooling tower is %lf pa or %lf inches of water col", tpd.metric, tpd.us);

    bhp.us = (cfm_fan.us * (tpd_s_inh20 + pd_v.us)) / (6356 * fann); //hp
    bhp.metric = bhp.us * 0.745699872; //kw
    printf("\n the required brake horse power of the fan is %lf hp or %lf kw", bhp.us, bhp.metric);

    return (bhp.metric);
}

I did have to make a small change to the data file, as well:

43.0 40.0 37.0 40.0 7.0 1.0 0.75
3.0 8.0 8.0 1000.0
3 4
273.15 273.16
-2948.997118 -2.1836674 -0.000150474 -0.0303738468 0.00042873 4.76955 25.83220018 -0.5
18.01534 28.9645 101325 0.000666
1.00568 1.84598 2500.84 8.31432
4 11 300 750

HTH.

Edited 11 Years Ago by Schol-R-LEA

Schol-R-LEA 1,446 Commie Mutant Traitor

11 Years Ago

OK, Let's trace that value backwards... fill_s[0] is ten to the power of xdd; xdd is the inverse of xd ; xd is the one minus pd_fill_n... where does the value of pd_fill_n come from?

I can only assume that pd_fill_n is a global, as it isn't declared in the code shown here. Is pd_fill_n initialized? if it isn't, then if your are lucky, the compiler is setting it to zero. More likely, it contains whatever garbage happens to be on the stack where the variable is stored. To fix fill_s, make sure that pd_fill_n has the correct value set.

Edited 11 Years Ago by Schol-R-LEA

Reply to this topic

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 1 · 2013-04-08T15:16:34+00:00

I am certain someone can, if you you show us what you've already done on the problem. If you haven't started working on it yet yourself, then I would recommend familiarizing yourself with the details of the algorithm before proceeding.

We are here to help you learn. This is not the same thing as helping you get a passing grade, though heopefully we can help you do that as well. We will help fix problems with your code, but we won't provide code for you.

sai.aleenus 0 Newbie Poster · Answer 2 · 2013-04-09T13:54:00+00:00

yeah..i have taken up the code for the algorithm but am facing some error problems,,
the code is as follows where i have taken out the bhpeqns func to facilitate for soace but its the major func to be optimised

#include <stdio.h>
#include <stdlib.h>
#include<conio.h>
#include<math.h>
#include<float.h>
#include<stdarg.h>
//#include<graphics.h>
#include<iostream.h>;


float bhpeqns(float x[vars],int n);

double hooke(int nvars,double startpt[vars],double endpt[vars],double rho,double epsilon,int itermax);

float bestsearch(double delta[vars],double point[vars],double prevbest,int nvars);



int vars,itermax;


float thw=43.0f,tcw=40.0f,twbi=37.0f,tdbi=40.0f,app,r,fand=7.0f,hubd=1.0f,fann=0.75f;
float aia,aih=3.0f,edisa,edisa_f,cell=8.0f,celw=8.0f,cela,twbo,tdbo,gpm=1000.0f,gpmpc,wl,spe_in,wl_f;
int ais = 3,celn = 4;
float to = 273.15, t9 = 273.16, svpa = -2948.997118, svpb = -2.1836674, svpc = -0.000150474,svpd = -0.0303738468;
float svpe = 0.00042873, svpf = 4.76955,svpg = 25.83220018;
float sphum1 = 18.01534, sphum2 = 28.9645, sphupo = 101325, sphuc = 0.000666;
float  spec2 = 1.00568, spec3 = 1.84598, spelo = 2500.84,spvolro = 8.31432;float row_fill_f;
int i,j,wlmin=4 , wlmax = 11, vmin = 300, vmax = 750;
int k;
// float wl_f,pdf1,pdf2,pdf3,pdf4,pdf5,pdf6;


#define vars (250)
#define rhoi (0.85)
#define epsimin (0.0001)
#define imax (1000)
#define n 2
#define nvars 2

int funevals = 0;


void main()
{
    float startpt[vars], endpt[vars];
    int nvars;int itermax;
    float rho, epsilon;
    int i, jj;

       /* starting guess for rosenbrock test function */
       nvars = 2;
       startpt[0] = 0.6892;
       startpt[1] = 12.0f;
       itermax = imax;
       rho = rhoi;
       epsilon = epsimin;
       jj = hooke(nvars, startpt, endpt, rho, epsilon, itermax);
       printf("\n\n\nHOOKE USED %d ITERATIONS, AND RETURNED\n", jj);
       for (i = 0; i < nvars; i++)
          {
             printf("x[%3d] = %15.7le \n", i, endpt[i]);
          }

}


float bhpeqns(float x[vars],int n)
{

    printf("i got called");


return (bhp_w);
}

float bestsearch(double delta[vars],double point[vars],double prevbest,int nvars)
{
       double z[vars],minf,ftmp;
       int i;

        minf = prevbest;
for(i = 0; i < nvars ; i++)
{
    z[i] = point[i];
}
    for(i=0 ; i < nvars ; i++)
    {
        z[i] = point[i]+delta[i] ;
        ftmp = bhpeqns(z , nvars);
        if(ftmp < minf)
        {
            minf = ftmp;
        }
        else
        {
            delta[i] = 0.0 - delta[i];
            z[i] = point[i]+delta[i] ;
            ftmp = bhpeqns(z , nvars);
            if(ftmp < minf)
            {
                minf = ftmp;
            }
            else
            {
                z[i] = point[i];
            }
        }
    }
for(i = 0; i < nvars ; i++)
{
     point[i] = z[i];
}

return(minf);
}

double hooke(int nvars,double startpt[vars],double endpt[vars],float rho,float epsilon,int itermax)
{
       double delta[vars];
       double newf, fbefore, steplength, tmp;
       double xbefore[vars], newx[vars];
       int i, j, keep;
       int iters, iadj;
       for (i = 0; i < nvars; i++)
       {
           newx[i] = xbefore[i] = startpt[i];
           delta[i] = fabs(startpt[i] * rho);
           if (delta[i] == 0.0)
           {
               delta[i] = rho;
           }
       }

       iadj = 0;
       steplength = rho;
       iters = 0;
       fbefore = bhpeqns(newx, nvars);
       newf = fbefore;
       while ((iters < itermax) && (steplength > epsilon))
        {
           iters++;
           iadj++;
           printf("\nAfter %5d funevals, f(x) =  %.4le at\n", funevals, fbefore);

           for (j = 0; j < nvars; j++)
           {
               printf("   x[%2d] = %.4le\n", j, xbefore[j]);

           }
              for (i = 0; i < nvars; i++)
           {
               newx[i] = xbefore[i];
           }

           newf = bestsearch(delta, newx, fbefore, nvars);
           keep = 1;
              while ((newf < fbefore) && (keep == 1))
             {
               iadj = 0;
               for (i = 0; i < nvars; i++)
                 {
                   if(newx[i] <= xbefore[i])
                    {
                     delta[i] = 0.0 - fabs(delta[i]);
                    }
                    else
                    {
                       delta[i] = fabs(delta[i]);
                    }
                   tmp = xbefore[i];
                   xbefore[i] = newx[i];
                   newx[i] = newx[i] + newx[i] - tmp;
               }

               fbefore = newf;
               newf = bestsearch(delta, newx, fbefore, nvars);

               if (newf >= fbefore)
                {
                   break;
                 }
               keep = 0;
               for (i = 0; i < nvars; i++)
                {
                   keep = 1;
                   if (fabs(newx[i] - xbefore[i]) > (0.5 * fabs(delta[i])))
                    {
                      break;
                     }
                   else
                     {
                       keep = 0;
                     }
               }
           }

             if ((steplength >= epsilon) && (newf >= fbefore))
             {
               steplength = steplength * rho;
               for (i = 0; i < nvars; i++)
               {
                   delta[i] *= rho;
               }
             }
       }
    for (i = 0; i < nvars; i++)
    {
          endpt[i] = xbefore[i];
    }
       return (iters);
}

the code is showing some lvaue erros and some others as well...can someone please help me with the error rectification...

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 3 · 2013-04-09T14:47:47+00:00

OK, before I address the issues of the functions themselves, I'd like to look at some technical issues in the code overall.

The main() function in both C and C++ should always return an int value; void main() is incorrect except in certain cases in C (and always in C++).
The <iostream.h> header is not only obsolete, it is specific to C++; if this is supposed to be a C program, it should not be included here.
The headers <conio.h> and <graphics.h> are specific to the Borland compilers; given that you don't actually use them anywhere, you can safely take them out.
Similarly, there is no standard <float.h> header, so this can also be removed.
Since you don't use variadic functions anywhere, <stdarg.h> can be removed as well.
You generally want to put all #define directives before any other code that references the constants. This relates to the next part, where you use a defined constant before its definition.
In the prototypes of the function bhpeqns(), you declare a parameter float x[vars], where vars has not yet been declared, never mind defined. In a declaration such as this, the array size must be a compile-time constant; if you don't know the size of the array, use float x[] instead. This is repeated in the declarations of hooke() and bestsearch() as well.
You have #defined the value of nvar as 2, but then you delcare an int nvar; in the main() function. Since the preprocessor replacements take place before the compiler handles declarations, this becomes int 2;, which is a syntax error.
You need to be more consistent in your indentation. The specific indent style you use is less important than being consistent in applying it.

I'll let you address these issues for now, while I look into the function code itself.

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 4 · 2013-04-09T19:12:26+00:00

At the risk of being presumptuous, I have taken the liberty of breaking the program into separate units, to make it easier to debug the separate parts. I have re-written the functions so that a) the function to be optimized is passed as an argument to the hooke() and bestsearch() functions, b) the vars value are passed to all three functions, and c) the hooke() and bestsearch() functions dynamically allocate their internal arrays, and free them when finished. While the code will compile and run, I do not know for certain if it is correct; you'll need to test it yourself for the correct results.

hooke-jeeves.h

#ifndef HOOKE_JEEVES_H
#define HOOKE_JEEVES_H 1

double hooke(double (* func)(double[], int, const int), const int vars, const int nvars, double startpt[], double endpt[], double rho, double epsilon, int itermax);

#endif

hooke-jeeves.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "hooke-jeeves.h"


double bestsearch(double (* func)(double[], int, const int), double delta[],double point[],double prevbest, const int vars, const int nvars);

int funevals = 0;


double hooke(double (* func)(double[], int, const int), const int vars, const int nvars, double startpt[], double endpt[], double rho, double epsilon, int itermax)
{
    double *delta;
    double newf, fbefore, steplength, tmp;
    double *xbefore, *newx;
    int i, j, keep;
    int iters, iadj;

    delta = calloc(vars, sizeof(double));

    if (delta == NULL)
    {
        printf("hooke: Memory allocation error variable delta\n");
        exit(-1);
    }

    xbefore = calloc(vars, sizeof(double));
    if (xbefore == NULL)
    {
        printf("hooke: Memory allocation error variable xbefore\n");
        free(delta);
        exit(-1);
    }

    newx = calloc(vars, sizeof(double));

    if (newx == NULL)
    {
        printf("hooke: Memory allocation error variable newx");
        free(delta);
        free(xbefore);
        exit(-1);
    }

    for (i = 0; i < nvars; i++)
    {
        newx[i] = xbefore[i] = startpt[i];
        delta[i] = fabs(startpt[i] * rho);
        if (delta[i] >= -epsilon && delta[i] <= epsilon)
        {
            delta[i] = rho;
        }
    }

    iadj = 0;
    steplength = rho;
    iters = 0;
    fbefore = func(newx, nvars, vars);
    newf = fbefore;
    while ((iters < itermax) && (steplength > epsilon))
    {
        iters++;
        iadj++;
        printf("\nAfter %5d funevals, f(x) =  %.4e at\n", funevals, fbefore);

        for (j = 0; j <  nvars; j++)
        {
            printf("   x[%2d] = %.4e\n", j, xbefore[j]);

        }
        for (i = 0; i < nvars; i++)
        {
            newx[i] = xbefore[i];
        }

        newf = bestsearch(func, delta, newx, fbefore, vars, nvars);
        keep = 1;
        while ((newf < fbefore) && (keep == 1))
        {
            iadj = 0;
            for (i = 0; i < nvars; i++)
            {
                if(newx[i] <= xbefore[i])
                {
                    delta[i] = 0.0 - fabs(delta[i]);
                }
                else
                {
                    delta[i] = fabs(delta[i]);
                }
                tmp = xbefore[i];
                xbefore[i] = newx[i];
                newx[i] = newx[i] + newx[i] - tmp;
            }

            fbefore = newf;
            newf = bestsearch(func, delta, newx, fbefore, vars, nvars);

            if (newf >= fbefore)
            {
                break;
            }
            keep = 0;
            for (i = 0; i < nvars; i++)
            {
                keep = 1;
                if (fabs(newx[i] - xbefore[i]) > (0.5 * fabs(delta[i])))
                {
                    break;
                }
                else
                {
                    keep = 0;
                }
            }
        }

        if ((steplength >= epsilon) && (newf >= fbefore))
        {
            steplength = steplength * rho;
            for (i = 0; i < nvars; i++)
            {
                delta[i] *= rho;
            }
        }
    }
    for (i = 0; i < nvars; i++)
    {
        endpt[i] = xbefore[i];
    }

    free(delta);
    free(xbefore);
    free(newx);
    return (iters);
}




double bestsearch(double (* func)(double[], int, const int), double delta[],double point[],double prevbest, const int vars, const int nvars)
{
    double *z, minf, ftmp;
    int i;

    z = calloc(sizeof(double), vars);

    if (z == NULL)
    {
        printf("bestsearch: Memory allocation error variable z");
        exit(-1);
    }

    minf = prevbest;
    for(i = 0; i < nvars; i++)
    {
        z[i] = point[i];
    }
    for(i=0 ; i < nvars; i++)
    {
        z[i] = point[i]+delta[i] ;
        ftmp = func(z , nvars, vars);
        funevals++;

        if(ftmp < minf)
        {
            minf = ftmp;
        }
        else
        {
            delta[i] = 0.0 - delta[i];
            z[i] = point[i] + delta[i] ;
            ftmp = func(z , nvars, vars);
            funevals++;

            if(ftmp < minf)
            {
                minf = ftmp;
            }
            else
            {
                z[i] = point[i];
            }
        }
    }
    for(i = 0; i < nvars; i++)
    {
        point[i] = z[i];
    }

    free(z);

    return(minf);
}

hooke-jeeves-test.c

#include <stdio.h>
#include "hooke-jeeves.h"


#define VARS (250)
#define RHO (0.85)
#define EPSILON (0.0001)
#define ITERMAX (1000)
#define NVARS (2)

double bhpeqns(double x[], int n, const int vars);
double banana(double x[], int n, const int vars);

int main()
{
    double startpt[VARS], endpt[VARS];
    int i, jj;

    /* starting guess for rosenbrock test function */
    startpt[0] = 0.6892;
    startpt[1] = 12.0f;
    jj = hooke(bhpeqns, VARS, NVARS, startpt, endpt, RHO, EPSILON, ITERMAX);
    printf("\n\n\nHOOKE USED %d ITERATIONS, AND RETURNED\n", jj);
    for (i = 0; i < NVARS; i++)
    {
        printf("x[%3d] = %15.7e \n", i, endpt[i]);
    }

    return 0;
}


double bhpeqns(double x[], int n, const int vars)
{

    printf("I got called ");


    return (1.0);
}


double banana(double x[], int n, const int vars)
{

   double a, b, c;

   a = x[0];
   b = x[1];
   c = 100.0 * (b - (a * a)) * (b - (a * a));
   return (c + ((1.0 - a) * (1.0 - a)));
}

I know you may not have encountered some of these constructs before, so let me know if you need further explanation.

sai.aleenus 0 Newbie Poster · Answer 5 · 2013-04-10T18:16:47+00:00

hey thanks for the suggestions and they were very help ful...can you also help me with the following code where the min of bhp is to be calculated while keeping the fh and fs which are the arguments passed in bhpeqns with in their ranges of fhi<=fhx[]<=fhl and fsi<=fsx[]<=fsl..
the code written by me is as follows where only the values in the immediate area of fh and fs are varying but not the entite range variation is being done...

#include <stdio.h>
#include <stdlib.h>
#include<conio.h>
#include<math.h>
#include<float.h>
#include<stdarg.h>
//#include<graphics.h>
#include<iostream.h>;



float bhpeqns(float fh,float fs);
float hooke(int max,float fh,float fs,float epsilon,float fhsl);



float thw=43.0f,tcw=40.0f,twbi=37.0f,tdbi=40.0f,app,r,fand=7.0f,hubd=1.0f,fann=0.75f;
float aia,aih=3.0f,edisa,edisa_f,cell=8.0f,celw=8.0f,cela,twbo,tdbo,gpm=1000.0f,gpmpc,wl,spe_in,wl_f;
int ais = 3,celn = 4;
float to = 273.15, t9 = 273.16, svpa = -2948.997118, svpb = -2.1836674, svpc = -0.000150474,svpd = -0.0303738468;
float svpe = 0.00042873, svpf = 4.76955,svpg = 25.83220018;
float sphum1 = 18.01534, sphum2 = 28.9645, sphupo = 101325, sphuc = 0.000666;
float  spec2 = 1.00568, spec3 = 1.84598, spelo = 2500.84,spvolro = 8.31432;float row_fill_f;
int i,j,wlmin=4 , wlmax = 11, vmin = 300, vmax = 750;
int k;
// float wl_f,pdf1,pdf2,pdf3,pdf4,pdf5,pdf6;



#define fhsl_ini 0.85
#define epsimin 0.0001
#define imax 100






void main()
{
  float fh,fs,fhsl,epsilon;
  int ima,k;
  printf("\n enter a value for fill height(note in b/w 0.6096 and 1.8288) : \t");
  scanf("%f",&fh);
  printf("\n enter a value for flute spacing(note in b/w 12 and 19) : \t");
  scanf("%f",&fs);
  fh = fh;
  fs = fs;
  fhsl = fhsl_ini;
  ima = imax;
  epsilon = epsimin;
  k = hooke(ima,fh,fs,epsilon,fhsl);
  printf("\n the iterations done are %d",k);
getch();


}


float bhpeqns(float fh, float fs)
{

    printf("i got called");

    float fh_f,fs_f;
    float svp_in,svph,svpi,svpj,svpk, svpl, svpm, svpn,too;//sat vap pressure
    float sphu_in,sphua,sphub,sphud,hin_f;//sp humidity
    float too1,too2,too3,too4,twbi1,twbi2,twbi3,twbi4,svp1,svp2,svp3,svp4,t1,t2,t3,t4;
    float sphu1,sphu2,sphu3,sphu4,spe1,spe2,spe3,spe4,spe;
    float h1,h2,h3,h4,h11,h22,h33,h44;
    float speat1,speat2,speat3,speat4,ra,h1a,h1b,h2a,h2b,h3a,h3b,h4a,h4b;
    float kav1,kav2,lg;
    float kav2a,kav2b,kav2c,kav2d,edisa_f;
    float fvel,fvel_f,l_f,l_s,g_f,g_s,spvol_in_s,spvol_in_f,row_in_s,row_in_f;
    float too_o,twbo,tdbo,svp_out,sphu_out,spvol_out_s,spvol_out_f,hout_s,hout_f,tdbo_f,row_out_s,row_out_f;
    float vel_in_s,vel_in_f,vel_out_s,vel_out_f,vel_fill_s,vel_fill_f,row_fill_s,row_fill_f,spvol_fill_s,spvol_fill_f;
    float wl_f,pdf1,pdf2,pdf3,pdf4,pdf5,pdf6,pd_il_inh20,pd_il_pa,cfm_fan_s,cfm_fan_f,pd_de_inh20,pd_de_pa,pd_fill_inh20,pd_fill_pa;
    float tpd_s_inh20,tpd_s_pa,pd_v_f,pd_v_s,bhp,bhp_w,tpd_f,tpd_s;
    float sphure,sphuref,spvolref,spvolre,predp;
    float pdvf2,pdvf3,bhhp1,bhhp2;



            r = thw - tcw ;
            app = tcw - twbi ;
            cela = cell*celw ;
            aia = ais*celw*aih;
            edisa = (0.785*((fand*fand)-(hubd*hubd)));
            edisa_f = edisa*3.2808*3.2808;
            printf("\n edisa %f",edisa_f);
            gpmpc = gpm/celn ;
            wl = gpmpc/cela;
            wl_f = (gpmpc*4.402867)/(cela*3.2808*3.2808);// fpm


           too = twbi+to;
            svph = svpa/too;//
            svpi = svpb*log(too);
            svpj = svpd*(twbi-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg ;
            svp_in = pow(10,svpn);                //sat vap pressure calculated
            printf("\n sat vap p inlet is %f Pa",svp_in);

            sphua = sphum1/sphum2;//
            sphub = sphupo*sphuc*(tdbi-twbi);
            sphud = svp_in - sphub;
            sphu_in = (sphua*sphud)/(sphupo - sphud); // sp humidity is calculated
            printf("\n sp humidity inlet is %f kg/kg of dry air",sphu_in);

            spe_in = (spec2*tdbi)+(sphu_in*(spelo+(spec3*tdbi)));//sp enthalpy is calculated
            printf("\n inlet enthalpy is(kj/kg) %f",spe_in);
            spe = spe_in;

            hin_f = ((spe_in*1.8)/4.1858)+7.686 ;
            printf("\n inlet enthalpy is(btu/lb) %f",hin_f);

            lg = 1.0f;
            do
            {

            t1 = tcw + (0.1*r);
            t2 = tcw + (0.4*r);
            t3 = thw - (0.4*r);
            t4 = thw - (0.1*r);

            twbi1 = twbi + (0.1*r);
            twbi2 = twbi + (0.4*r);
            twbi3 = twbi + (0.6*r);
            twbi4 = twbi + (0.9*r);

            too1 = twbi1+to;
            svph = svpa/too1;//
            svpi = svpb*log(too1);//
            svpj = svpd*(twbi1-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too1));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg ;
            svp1 = pow(10,svpn);                //sat vap pressure calculated
         //   printf("\n sat vap p1 is %f",svp1);

                                  //too1,twbi1,svp1,sphu1,spe1
            sphua = sphum1/sphum2;//
            sphub = sphupo*sphuc*(t1-twbi1);
            sphud = svp1 - sphub;
            sphu1 = (sphua*sphud)/(sphupo - sphud); // sp humidity is calculated
          //  printf("\n sp humidity1 is %f",sphu1);


            spe1 = (spec2*t1)+(sphu1*(spelo+(spec3*t1)));//sp enthalpy is calculated
          //  printf("\n enthalpy1 is %f",spe1);

            too2 = twbi2+to;
            svph = svpa/too2;//
            svpi = svpb*log(too2);//
            svpj = svpd*(twbi2-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too2));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg ;
            svp2 = pow(10,svpn);                //sat vap pressure calculated
          //  printf("\n sat vap p2 is %f",svp2);

                                  //
            sphua = sphum1/sphum2;//
            sphub = sphupo*sphuc*(t2-twbi2);
            sphud = svp2 - sphub;
            sphu2 = (sphua*sphud)/(sphupo - sphud); // sp humidity is calculated
          //  printf("\n sp humidity2 is %f",sphu2);


            spe2 = (spec2*t2)+(sphu2*(spelo+(spec3*t2)));//sp enthalpy is calculated
           // printf("\n enthalpy2 is %f",spe2);



            too3 = twbi3+to;
            svph = svpa/too3;//
            svpi = svpb*log(too3);//
            svpj = svpd*(twbi3-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too3));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg ;
            svp3 = pow(10,svpn);                //sat vap pressure calculated
          //  printf("\n sat vap p3 is %f",svp3);

                                  //too1,twbi1,svp1,sphu1,spe1
            sphua = sphum1/sphum2;//
            sphub = sphupo*sphuc*(t3-twbi3);
            sphud = svp3 - sphub;
            sphu3 = (sphua*sphud)/(sphupo - sphud); // sp humidity is calculated
           // printf("\n sp humidity3 is %f",sphu3);

            spe3 = (spec2*t3)+(sphu3*(spelo+(spec3*t3)));//sp enthalpy is calculated
          //  printf("\n enthalpy3 is %f",spe3);


            too4 = twbi4+to;
            svph = svpa/too4;//
            svpi = svpb*log(too4);//
            svpj = svpd*(twbi4-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too4));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg ;
            svp4 = pow(10,svpn);                //sat vap pressure calculated
            //printf("\n sat vap p4 is %f",svp4);

                                  //too1,twbi1,svp1,sphu1,spe1
            sphua = sphum1/sphum2;//
            sphub = sphupo*sphuc*(t4-twbi4);
            sphud = svp4 - sphub;
            sphu4 = (sphua*sphud)/(sphupo - sphud); // sp humidity is calculated
           // printf("\n sp humidity4 is %f",sphu4);


            spe4 = (spec2*t4)+(sphu4*(spelo+(spec3*t4)));//sp enthalpy is calculated
           // printf("\n enthalpy4 is %f",spe4);

            speat1 = spe + (0.1*lg*r);
            speat2 = spe + (0.4*lg*r);
            speat3 = spe + (0.6*lg*r);
            speat4 = spe + (0.9*lg*r);

            h1a = spe1-speat1;
            h1b = ((h1a*1.8)/4.1858);
            h1 = h1b + 7.686 ;
            h2a = spe2-speat2;
            h2b = ((h2a*1.8)/4.1858);
            h2 = h2b + 7.686 ;
            h3a = spe3-speat3;
            h3b = ((h3a*1.8)/4.1858);
            h3 = h3b + 7.686 ;
            h4a = spe4-speat4;
            h4b = ((h4a*1.8)/4.1858);
            h4 = h4b + 7.686 ;//btu/lb
            h11 = 1/h1;
            h22 = 1/h2;
            h33 = 1/h3;
            h44 = 1/h4;

            ra = (r*1.8)+32;//in deg f

            kav1 = ra*0.25*(h11+h22+h33+h44);//in fps ---- demand

 //    fvel = (wl*8.33)/(lg*0.07);//in m/s
            fvel_f = (wl*8.33)/(lg*0.07*0.227125);//in fpm conversion factor 0.227125 used

            fh_f = fh/0.3048;
            fs_f = fs*0.00328084;//in ft

            kav2a = 0.06717927* pow(lg,-0.546673);
            kav2b = pow(fvel_f,.25158684);
            kav2c = pow(fh_f,0.6369974);
            kav2d = pow (fs_f,-0.429);
            kav2 = kav2a*kav2b*kav2c*kav2d; //supply ntu

           /* for(lg=1; kav1<=kav2;lg+0.1)
            {
            printf("\n the l/g ratio required for the performance is %f \n",lg);
            printf("\n the KAV/L  for the performance is %f \n",kav2);
            }*/
            /*lg = 1;*/


            //printf("\n the l/g ratio required for the performance is %f \n",lg);
            lg = lg+0.001;
            continue;
               /* if(lg2<lg)
            {
                lg = lg2 + 0.001;

            }
            else{lg=lg+0.2;}*/

            }while(kav1<=kav2 && lg<=3 && lg>=1);
           /* while(lg=1 && lg<=5)
            {
            if(kav2>=kav1)
            break;
            else
            lg = lg+0.2;
            printf("\n the l/g ratio required for the performance is %f \n",lg);
            }*/

               printf("\n the l/g ratio required for the performance is %f \n",lg);
               kav2 = kav2a*kav2b*kav2c*kav2d;
               printf("\n the KAV/L_2  for the performance is %f \n",kav2);
               kav1 = ra*0.25*(h11+h22+h33+h44);
               printf("\n the KAV/L_1  for the performance is %f \n",kav1);

 // fvel = (wl*8.33)/(lg*0.07);//in m/s
               fvel_f = (wl*8.33)/(lg*0.07*0.00508);
               //printf("\n the velocity at the fill is %f m/s or %f fpm \n",fvel,fvel_f);

               l_f = (gpmpc *8.33)/0.227125;//lb/min
               l_s = (gpmpc * 1000)/3600;// kg/sec
               l_f = l_s/0.00755987283;
               printf("\n the mass flow rate of water is %f kg/s or %f lb/min \n",l_s,l_f);

               g_f = l_f/lg;
               g_s = l_s/lg;
               printf("\n the mass flow rate of air is %f kg/s or %f lb/min \n",g_s,g_f);



               printf("\n AT COOLING TOWER INLET\n");

               printf("\n sat vap p inlet is %f Pa",svp_in);
               printf("\n sp humidity inlet is %f kg/kg of dry air",sphu_in);
               spvol_in_s = (spvolro*(tdbi+to)*1000)/(sphum2*(sphupo-svp_in+(sphupo*sphuc*(tdbi-twbi))));//cum/kg
               spvol_in_f = spvol_in_s / 0.0624279606;//cuft/lb
               printf("\n the inlet specific volume is %f cum/kg or %f cu ft/lb ",spvol_in_s,spvol_in_f);//

               row_in_s =(1+ sphu_in)/(spvol_in_s);
               row_in_f = row_in_s/16.0184634;
               printf("\n the inlet density is %f kg/cum or %f lb/cu ft",row_in_s,row_in_f);

               hout_s = spe + (lg*r);
               hout_f = ((hout_s*1.8)/4.1858)+7.686 ;
               printf("\n the inlet enthalpy is %f kj/kg or %f btu/lb",spe_in,hin_f);
               printf("\n the exit enthalpy is %f kj/kg or %f btu/lb",hout_s,hout_f);

               tdbo = hout_s/1.006;
               printf("\n tdbo is %f",tdbo);
           //    tdbo_f = hout_f/0.240;//deg f
               twbo = tdbo;

            too_o = twbo+to;
            svph = svpa/too_o;//
            svpi = svpb*log(too_o);                     //
            svpj = svpd*(twbo-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too_o));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg-0.5 ;
            svp_out = pow(10,svpn);                //sat vap pressure calculated
            printf("\n sat vap p exit is %f Pa   svpn    %f",svp_out,svpn);

            sphua = sphum1/sphum2;//
            sphub = 0;
            sphud = svp_out ;
            sphure = sphud - sphupo ;
            sphuref = abs(sphure);

            sphu_out = (sphua*sphud)/(sphuref);
//          sphuout = (sphua*sphud)/(sphupo - sphud);
            printf("\n sp humidity exit is %f kg/kg of satureated air",sphu_out);


            printf("\n AT COOLING TOWER EXIT\n");

               spvolre = svp_out - sphupo;
            spvolref = abs(spvolre);
               spvol_out_s = (spvolro*(tdbo+to)*1000)/(sphum2*(spvolref));//cum/kg
               spvol_out_f = spvol_out_s / 0.0624279606;//cuft/lb
               printf("\n the exit specific volume is %f cum/kg or %f cu ft/lb ",spvol_out_s,spvol_out_f);

               row_out_s =(1+ sphu_out)/(spvol_out_s);

               row_out_f = row_out_s/16.0184634;
               printf("\n the exit density is %f kg/cum or %f lb/cu ft",row_out_s,row_out_f);

               vel_in_s = (g_s * spvol_in_s)/aia ;//m/s
               vel_in_f = vel_in_s/0.00508;//fpm
               vel_out_s = (g_s*spvol_out_s)/cela ;
               vel_out_f =  vel_out_s/0.00508;

               row_fill_s = (row_in_s + row_out_s)/2 ;
               row_fill_f = (row_in_f + row_out_f)/2 ;
               spvol_fill_s = (spvol_in_s+spvol_out_s)/2;
               spvol_fill_f = (spvol_in_f+spvol_out_f)/2;   //

               vel_fill_s = (g_s*spvol_fill_s)/cela ;
               vel_fill_f = vel_fill_s/0.00508;

               printf("\n the specific volume at fill is %f cum/kg or %f cu ft/lb ",spvol_fill_s,spvol_fill_f);
               printf("\n the density at fill is %f kg/cum or %f lb/cu ft",row_fill_s,row_fill_f);
               printf("\n the velocity of air at inlet is %f m/s or %f fpm ",vel_in_s,vel_in_f);
               printf("\n the velocity of air at fill is %f m/s or %f fpm ",vel_fill_s,vel_fill_f);
               printf("\n the velocity of air at exit is %f m/s or %f fpm ",vel_out_s,vel_out_f);

               cfm_fan_s = g_s * spvol_out_s;//cum/s
               cfm_fan_f = g_f * spvol_out_f;
               printf("\n the volume flow rate of air at fan CFAN is %f cum/s or %f cu ft/min \n",cfm_fan_s,cfm_fan_s);



               pd_il_inh20 = (3*((2*pow(10,-7)*vel_in_f*vel_in_f)+(3*pow(10,-6)*vel_in_f))*(row_in_f/0.07));
        //       pd_il_pa = pd_il_inh20*249;//pascal
               printf("\n pd il %f in ",pd_il_inh20);

               pd_de_inh20 = (5* ((3*pow(10,-7)*vel_out_f*vel_out_f)-(9*pow(10,-6)*vel_out_f) + 0.003 )*(row_out_f/0.07));
        //       pd_de_pa = pd_de_inh20*249;//pascal
               printf("\n pd de %f in ",pd_de_inh20);



               pdf1 = (fh_f/fs_f);
               pdf2 =pow(pdf1,0.963077);
               pdf3 = ((7.44385*pow(10,-6))+((-4.078*pow(10,-6))*pow(fh_f,0.176686)));
               pdf4 = (4.897*(pow(10,-4))*(pow(vel_fill_f,2.3386)));
               pdf5 = (0.0428*wl_f);
               pdf6 = (pow(vel_fill_f,1.33558));

               pd_fill_inh20 = (pdf2*pdf3*(pdf4+(pdf5*pdf6)))*(row_fill_f/0.07);
               printf("\n pd fill %f in ",pd_fill_inh20);
         //      pd_fill_pa = pd_fill_inh20*249;

               tpd_s_inh20 = pd_il_inh20 + pd_de_inh20 + pd_fill_inh20;
               printf("\n hs %f in ",tpd_s_inh20);
               tpd_s_pa = tpd_s_inh20*249;
               pd_v_f = (row_out_f/1202312)*(cfm_fan_f/edisa_f)*(cfm_fan_f/edisa_f);
               printf("\n hv %f in",pd_v_f);
               pd_v_s = pd_v_f * 249 ;
               tpd_f = pd_v_f + tpd_s_inh20;
               tpd_s = tpd_s_pa + pd_v_s;
               printf("\n the overall pressure drop across the cooling tower is %f pa or %f inches of water col",tpd_s,tpd_f);

               bhp = ((cfm_fan_f*(tpd_s_inh20+pd_v_f))/(6356*fann));//hp
               bhp_w = bhp*0.745699872;//kw
               printf("\n the required brake horse power of the fan is %f Bhp or %f kw",bhp,bhp_w);

return (bhp_w);
}



float hooke(int max,float fh,float fs,float epsilon,float fhsl)
{

int iadj,iters,i,j,xs1,ys1,xs2,ys2,islfs = 1,fssl = 1;
float islfh,fhi,fhf,fsi,fsf,bhpi;

float a,b,m;
float bhp[250],fh1ta[250],fs1ta[250],bhp1ta[250],fh1pa,fs1pa,bhp1pa;
float fh1tn[250],fs1tn[250],bhp1tn[250],fh1pn,fs1pn,bhp1pn,fh1[250],fs1[250],bhp1[250],fh1_p,fs1_p,bhp1_p,fh1_m,fs1_m,bhp1_m;
float fh2ta[250],fs2ta[250],bhp2ta[250],fh2tn[250],fs2tn[250],bhp2tn[250];
float fh2pa,fs2pa,bhp2pa,fh2pn,fs2pn,bhp2pn,fh2_m,fs2_m,bhp2_m,fh2[250],fs2[250],bhp2[250];
float fhsp[250],fssp[250],lam1[250],lam2[250],lam[250],delta,fh3[250],fs3[250],bhp3[250],fh3_m,fs3_m,bhp3_m;


islfh = 0.001;
xs1 = 1;
ys1 = 0;
xs2 = 0;
ys2 = 1;
iadj = 0;
iters = 0;
fhi = 0.6096;
fhf = 1.8288;
fsi = 12.0f;
fsf = 19.0f;
a = (7.44385*pow(10,-6));
b = (-4.078*pow(10,-6));
m = 0.176686;


bhpi = bhpeqns(fhi , fsi);
while( (fh >= fhi) && (fh<= fhf) && (fs >= fsi) && (fs <= fsf) && (fhsl >= epsilon))
{
    iadj++;
    fhi = fh ; fsi = fs;
    for(i = 0 ; (i < max) && (bhp[i+1] <= bhp[i]) ; i++)
    {
        iters++;
        fh1ta[i] = fhi + (xs1 * islfh);
        fs1ta[i] = fsi + (ys1 * islfs);
        fh1pa = fh1ta[i] ;      fs1pa = fs1ta[i] ;
        bhp1pa = bhpeqns(fh1pa,fs1pa);
        bhp1ta[i] = bhp1pa ;

        fh1tn[i] = fhi - (xs1 * islfh);
        fs1tn[i] = fsi - (ys1 * islfs);
        fh1pn = fh1tn[i] ;      fs1pn = fs1tn[i] ;
        bhp1pn = bhpeqns(fh1pn,fs1pn);
        bhp1tn[i] = bhp1pn ;

        if((bhp1pa < bhpi) && (fh1[i] >= fhi) && (fh1[i] <= fhf) && ( fs1[i] >= fsi) && (fs1[i] <= fsf))
        {
            for(j = 0 ; (bhp1[j+1]<= bhp1[j])  && (fhsl >= islfh) ; j++)
            {
                fh1[j] = fhi + (xs1 * fhsl);
                fs1[j] = fsi + ( ys1 * fssl);
                fh1_p = fh1[j]  ;   fs1_p = fs1[j] ;
                bhp1_p = bhpeqns(fh1_p,fs1_p);
                bhp1[j] = bhp1_p;

                if(bhp1_p < bhpi)
                {
                   fh1[j] = fhi + (xs1 * fhsl);
                   fs1[j] = fsi + ( ys1 * fssl);
                   fh1_p = fh1[j] ;     fs1_p = fs1[j];
                }
                else
                {
                    fhsl = fhsl/2 ;
                    j++ ;
                }
            }
                fh1[i] = fh1_p;
                fs1[i] = fs1_p;

        }
        else if((bhp1pn < bhpi) && (fh1[i] >= fhi) && (fh1[i] <= fhf) && ( fs1[i] >= fsi) && (fs1[i] <= fsf))
        {
            for(j = 0 ; (bhp1[j+1]<= bhp1[j])  && (fhsl >= islfh) ; j++)
            {
                fh1[j] = fhi - (xs1 * fhsl);
                fs1[j] = fsi - ( ys1 * fssl);
                fh1_p = fh1[j] ;    fs1_p = fs1[j];
                bhp1_p = bhpeqns(fh1_p,fs1_p);
                bhp1[j] = bhp1_p;

                if(bhp1_p < bhpi)
                {
                   fh1[j] = fhi - (xs1 * fhsl);
                   fs1[j] = fsi - ( ys1 * fssl);
                   fh1_p = fh1[j] ;     fs1_p = fs1[j];
                }
                else
                {
                    fhsl = fhsl/2 ;
                    j++ ;
                }
            }
                fh1[i] = fh1_p;
                fs1[i] = fs1_p;

        }
        else
        {
            fh1[i] = fhi ;
            fs1[i] = fsi ;
        }

        fh1_m = fh1[i] ;
        fs1_m = fs1[i] ;
        bhp1_m = bhpeqns(fh1_m,fs1_m);
        bhp1[i] = bhp1_m ;
        fh = fh1[i];    fs = fs1[i];    bhp[i] = bhp1[i];
        printf("\n  i = %d  fh1   %f   fs1   %f    bhp1    %f",i,fh1[i],fs1[i],bhp1[i]);

        fh2ta[i] = fh1_m + (xs2 * islfh);
        fs2ta[i] = fs1_m + (ys2 * islfs);
        fh2pa = fh2ta[i];   fs2pa = fs2ta[i];
        bhp2pa = bhpeqns(fh2pa , fs2pa);
        bhp2ta[i] = bhp2pa;

        fh2tn[i] = fh1_m - (xs2 * islfh);
        fs2tn[i] = fs1_m - (ys2 * islfs);
        fh2pn = fh2tn[i];   fs2pn = fs2tn[i];
        bhp2pn = bhpeqns(fh2pn , fs2pn);
        bhp2tn[i] = bhp2pn;

        if((bhp2pa < bhp1_m)&& (fh2[i] >= fhi) && (fh2[i] <= fhf) && ( fs2[i] >= fsi) && (fs2[i] <= fsf))
        {
            fh2[i] = fh1_m + (xs2 * fhsl);
            fs2[i] = fs1_m + (ys2 * fssl);
        }
        else if( (bhp2pn < bhp1_m)  && (fh2[i] >= fhi) && (fh2[i] <= fhf) && ( fs2[i] >= fsi) && (fs2[i] <= fsf))
        {
            fh2[i] = fh1_m - (xs2 * fhsl);
            fs2[i] = fs1_m - (ys2 * fssl);
        }
        else
        {
            fh2[i] = fh1_m;
            fs2[i] = fs1_m;
        }

        fh2_m = fh2[i] ;
        fs2_m = fs2[i] ;
        bhp2_m = bhpeqns(fh2_m,fs2_m);
        bhp2[i] = bhp2_m = bhp[i];
        fh = fh2[i];    fs = fs2[i];    bhp[i] = bhp2[i];
        printf("\n  i = %d  fh2   %f   fs2   %f    bhp2    %f",i,fh2[i],fs2[i],bhp2[i]);

        fhsp[i] = fh2_m - fh1_m;
        fssp[i] = fs2_m - fs1_m;
        lam1[i] = (-1*a)/(b*fh2_m*m);
        lam2[i] = 1/(m-1);
        lam[i] = (pow(lam1[i],lam2[i]))-1 ;
        delta = lam[i];
        printf("\n the value of direction search i  %d   delta   %f",i,delta);

        fh3_m = fh2_m + (lam[i] * fhsp[i]);
        fh3[i] = fh3_m;
        fs3_m = fs2_m + (lam[i] * fssp[i]);
        fs3[i] = fs3_m;
        bhp3_m = bhpeqns(fh3_m,fs3_m);
        bhp3[i] = bhp3_m;
        bhp[i] = bhp3[i];   fh = fh3[i];    fs = fs3[i];
        printf("\n  i = %d  fh3   %f   fs3   %f    bhp3    %f",i,fh3[i],fs3[i],bhp3[i]);
        while((fh3[i] >= fhi) && (fh3[i] <= fhf) && ( fs3[i] >= fsi) && (fs3[i] <= fsf))
        {
        if (bhp3_m < bhp2_m)
        {
            fh3[i] = fh3[i];
            fs3[i] = fs3[i];
            bhp3[i] = bhp3[i];
            bhp3_m = bhp3[i];
            printf("\n the optimized value is bhp   %f   fh %f    fs    %f ",bhp3[i],fh3[i],fs3[i]);
            fhi = fh3_m;
            fh = fhi;
            fsi = fs3_m;
            fs = fsi;
            i++;
        }
        else
        {
            fh3[i] = fh2[i];
            fs3[i] = fs2[i];
            bhp3[i] = bhp2[i];
            bhp3_m = bhp3[i];
            printf("\n the optimized value is bhp   %f   fh %f    fs    %f ",bhp3_m,fh3[i],fs3[i]);
            fhi = fh2_m;
            fh = fhi;
            fsi = fs2_m;
            fs = fsi;
            i++;
        }
        }
        i++;
        printf("\n the number of times the iterations are run is %d",iters);
    }
    printf("\n the number of times the while loop is run is %d",iadj);

}
    for(i = 0 ; (i < imax) && (bhp[i+1] <= bhp[i]) ; i++)
    {

        printf("\n  i = %d  fh1   %f   fs1   %f    bhp1    %f",i,fh1[i],fs1[i],bhp1[i]);

        printf("\n  i = %d  fh2   %f   fs2   %f    bhp2    %f",i,fh2[i],fs2[i],bhp2[i]);
        printf("\n  i = %d  fh3   %f   fs3   %f    bhp3    %f",i,fh3[i],fs3[i],bhp3[i]);
    }
    printf("\n the optimized value is bhp   %f   fh %f    fs    %f ",bhp3_m,fh3[i],fs3[i]);
getch();
return (iters);
}

the major aim here is to vary the fh and fs values in their corresponding ranges to obtain the min of bhp and note that fs is to remain an integer at any cost...

can anyone resolve the issue for me as the ans is supposed to be fhf and fsi while i am getting very incorrect result..also can anyone suggest how i can represent the bhp variation with changes in fh and fs using graphs in c

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 6 · 2013-04-10T18:59:35+00:00

Ye gods and little fishes! How many different functions is bhpeqns() computing? This is insane! You must - absolutely must - break this down into separate functions if you stand any chance of understanding and debugging this mess. Trying to interleave multiple functions this way is courting disaster.

What you did with the Hookes-Jeeves function is particularly worrisome, since you took a fairly simple, modular function and made it into a morass of complication. This function should be something that you could apply to several different functions without even knowing the details of the functions it is applied to.

hey thanks for the suggestions and they were very help ful..

I can't see how helpful it was, given that you outright ignored the majority of my advice, and generally misinterpreted the rest.

sai.aleenus 0 Newbie Poster · Answer 7 · 2013-04-11T12:06:45+00:00

i have just tried this prrevious program and am reverting back to my first program code...
i have followed your guidelines and have corrected the errors and am finally getting the value of optimum bhp for the given arguments of fh and fs but values of fh and fs are skyrocketing.can you please help me by suggesting how i can constrain the x[0] and x[1] values for the first program of this discussion thread...
i want to implement 0.6096<x[0]<1.8288 and 12<x[0]<19..please help me with this

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 8 · 2013-04-11T13:14:59+00:00

Could you please post both your current code, and the problematic results? It would help if we were all looking at the same thing here.

sai.aleenus 0 Newbie Poster · Answer 9 · 2013-04-11T14:47:25+00:00

#include <stdio.h>
#include <stdlib.h>
#include<conio.h>
#include<math.h>
#include<float.h>
#include<stdarg.h>
//#include<graphics.h>
#include<iostream.h>;

#define vars 250
#define rhoi 0.85
#define epsimin (0.0001)
#define imax (1000)


float bhpeqns(float x[vars],int n);
float hooke(int nvars,float startpt[vars],float endpt[vars],float rho,float epsilon,int itermax);

float bestsearch(float delta[vars],float point[vars],float prevbest,int nvars);





float thw=43.0f,tcw=40.0f,twbi=37.0f,tdbi=40.0f,app,r,fand=7.0f,hubd=1.0f,fann=0.75f;
float aia,aih=3.0f,edisa,edisa_f,cell=8.0f,celw=8.0f,cela,twbo,tdbo,gpm=1000.0f,gpmpc,wl,spe_in,wl_f;
int ais = 3,celn = 4;
float to = 273.15, t9 = 273.16, svpa = -2948.997118, svpb = -2.1836674, svpc = -0.000150474,svpd = -0.0303738468;
float svpe = 0.00042873, svpf = 4.76955,svpg = 25.83220018;
float sphum1 = 18.01534, sphum2 = 28.9645, sphupo = 101325, sphuc = 0.000666;
float  spec2 = 1.00568, spec3 = 1.84598, spelo = 2500.84,spvolro = 8.31432;float row_fill_f;
int i,j,wlmin=4 , wlmax = 11, vmin = 300, vmax = 750;
int k;
// float wl_f,pdf1,pdf2,pdf3,pdf4,pdf5,pdf6;


//#define vars 250



int funevals = 0;

int main()
{
    float startpt[vars], endpt[vars];
    int nvars;int itermax;
    float rho, epsilon;
    int i, jj;

       /* starting guess for rosenbrock test function */
       nvars = 2;
       startpt[0] = 0.6892;
       startpt[1] = 12.0f;
       itermax = imax;
       rho = rhoi;
       epsilon = epsimin;
       jj = hooke(nvars, startpt, endpt, rho, epsilon, itermax);
       printf("\n\n\nHOOKE USED %d ITERATIONS, AND RETURNED\n", jj);
       for (i = 0; i < nvars; i++)
          {
             printf("x[%3d] = %15.7le \n", i, endpt[i]);
          }
      getch();
      return(i);

}


float bhpeqns(float x[vars], int n)
{

    printf("i got called");

    float fh_f,fs_f,fh,fs;
    float svp_in,svph,svpi,svpj,svpk, svpl, svpm, svpn,too;//sat vap pressure
    float sphu_in,sphua,sphub,sphud,hin_f;//sp humidity
    float too1,too2,too3,too4,twbi1,twbi2,twbi3,twbi4,svp1,svp2,svp3,svp4,t1,t2,t3,t4;
    float sphu1,sphu2,sphu3,sphu4,spe1,spe2,spe3,spe4,spe;
    float h1,h2,h3,h4,h11,h22,h33,h44;
    float speat1,speat2,speat3,speat4,ra,h1a,h1b,h2a,h2b,h3a,h3b,h4a,h4b;
    float kav1,kav2,lg;
    float kav2a,kav2b,kav2c,kav2d,edisa_f;
    float fvel,fvel_f,l_f,l_s,g_f,g_s,spvol_in_s,spvol_in_f,row_in_s,row_in_f;
    float too_o,twbo,tdbo,svp_out,sphu_out,spvol_out_s,spvol_out_f,hout_s,hout_f,tdbo_f,row_out_s,row_out_f;
    float vel_in_s,vel_in_f,vel_out_s,vel_out_f,vel_fill_s,vel_fill_f,row_fill_s,row_fill_f,spvol_fill_s,spvol_fill_f;
    float wl_f,pdf1,pdf2,pdf3,pdf4,pdf5,pdf6,pd_il_inh20,pd_il_pa,cfm_fan_s,cfm_fan_f,pd_de_inh20,pd_de_pa,pd_fill_inh20,pd_fill_pa;
    float tpd_s_inh20,tpd_s_pa,pd_v_f,pd_v_s,bhp,bhp_w,tpd_f,tpd_s;
    float sphure,sphuref,spvolref,spvolre,predp;
    float pdvf2,pdvf3,bhhp1,bhhp2;

    funevals++;
            printf("\n %d",n);
            r = thw - tcw ;
            app = tcw - twbi ;
            cela = cell*celw ;
            aia = ais*celw*aih;
            edisa = (0.785*((fand*fand)-(hubd*hubd)));
            edisa_f = edisa*3.2808*3.2808;
            printf("\n edisa %f",edisa_f);
            gpmpc = gpm/celn ;
            wl = gpmpc/cela;
            wl_f = (gpmpc*4.402867)/(cela*3.2808*3.2808);// fpm
        fh = x[0];
        fs = x[1];

           too = twbi+to;
            svph = svpa/too;//
            svpi = svpb*log(too);
            svpj = svpd*(twbi-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg ;
            svp_in = pow(10,svpn);                //sat vap pressure calculated
            printf("\n sat vap p inlet is %f Pa",svp_in);

            sphua = sphum1/sphum2;//
            sphub = sphupo*sphuc*(tdbi-twbi);
            sphud = svp_in - sphub;
            sphu_in = (sphua*sphud)/(sphupo - sphud); // sp humidity is calculated
            printf("\n sp humidity inlet is %f kg/kg of dry air",sphu_in);

            spe_in = (spec2*tdbi)+(sphu_in*(spelo+(spec3*tdbi)));//sp enthalpy is calculated
            printf("\n inlet enthalpy is(kj/kg) %f",spe_in);
            spe = spe_in;

            hin_f = ((spe_in*1.8)/4.1858)+7.686 ;
            printf("\n inlet enthalpy is(btu/lb) %f",hin_f);

            lg = 1.0f;
            do
            {

            t1 = tcw + (0.1*r);
            t2 = tcw + (0.4*r);
            t3 = thw - (0.4*r);
            t4 = thw - (0.1*r);

            twbi1 = twbi + (0.1*r);
            twbi2 = twbi + (0.4*r);
            twbi3 = twbi + (0.6*r);
            twbi4 = twbi + (0.9*r);

            too1 = twbi1+to;
            svph = svpa/too1;//
            svpi = svpb*log(too1);//
            svpj = svpd*(twbi1-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too1));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg ;
            svp1 = pow(10,svpn);                //sat vap pressure calculated
         //   printf("\n sat vap p1 is %f",svp1);

                                  //too1,twbi1,svp1,sphu1,spe1
            sphua = sphum1/sphum2;//
            sphub = sphupo*sphuc*(t1-twbi1);
            sphud = svp1 - sphub;
            sphu1 = (sphua*sphud)/(sphupo - sphud); // sp humidity is calculated
          //  printf("\n sp humidity1 is %f",sphu1);


            spe1 = (spec2*t1)+(sphu1*(spelo+(spec3*t1)));//sp enthalpy is calculated
          //  printf("\n enthalpy1 is %f",spe1);

            too2 = twbi2+to;
            svph = svpa/too2;//
            svpi = svpb*log(too2);//
            svpj = svpd*(twbi2-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too2));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg ;
            svp2 = pow(10,svpn);                //sat vap pressure calculated
          //  printf("\n sat vap p2 is %f",svp2);

                                  //
            sphua = sphum1/sphum2;//
            sphub = sphupo*sphuc*(t2-twbi2);
            sphud = svp2 - sphub;
            sphu2 = (sphua*sphud)/(sphupo - sphud); // sp humidity is calculated
          //  printf("\n sp humidity2 is %f",sphu2);


            spe2 = (spec2*t2)+(sphu2*(spelo+(spec3*t2)));//sp enthalpy is calculated
           // printf("\n enthalpy2 is %f",spe2);



            too3 = twbi3+to;
            svph = svpa/too3;//
            svpi = svpb*log(too3);//
            svpj = svpd*(twbi3-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too3));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg ;
            svp3 = pow(10,svpn);                //sat vap pressure calculated
          //  printf("\n sat vap p3 is %f",svp3);

                                  //too1,twbi1,svp1,sphu1,spe1
            sphua = sphum1/sphum2;//
            sphub = sphupo*sphuc*(t3-twbi3);
            sphud = svp3 - sphub;
            sphu3 = (sphua*sphud)/(sphupo - sphud); // sp humidity is calculated
           // printf("\n sp humidity3 is %f",sphu3);

            spe3 = (spec2*t3)+(sphu3*(spelo+(spec3*t3)));//sp enthalpy is calculated
          //  printf("\n enthalpy3 is %f",spe3);


            too4 = twbi4+to;
            svph = svpa/too4;//
            svpi = svpb*log(too4);//
            svpj = svpd*(twbi4-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too4));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg ;
            svp4 = pow(10,svpn);                //sat vap pressure calculated
            //printf("\n sat vap p4 is %f",svp4);

                                  //too1,twbi1,svp1,sphu1,spe1
            sphua = sphum1/sphum2;//
            sphub = sphupo*sphuc*(t4-twbi4);
            sphud = svp4 - sphub;
            sphu4 = (sphua*sphud)/(sphupo - sphud); // sp humidity is calculated
           // printf("\n sp humidity4 is %f",sphu4);


            spe4 = (spec2*t4)+(sphu4*(spelo+(spec3*t4)));//sp enthalpy is calculated
           // printf("\n enthalpy4 is %f",spe4);

            speat1 = spe + (0.1*lg*r);
            speat2 = spe + (0.4*lg*r);
            speat3 = spe + (0.6*lg*r);
            speat4 = spe + (0.9*lg*r);

            h1a = spe1-speat1;
            h1b = ((h1a*1.8)/4.1858);
            h1 = h1b + 7.686 ;
            h2a = spe2-speat2;
            h2b = ((h2a*1.8)/4.1858);
            h2 = h2b + 7.686 ;
            h3a = spe3-speat3;
            h3b = ((h3a*1.8)/4.1858);
            h3 = h3b + 7.686 ;
            h4a = spe4-speat4;
            h4b = ((h4a*1.8)/4.1858);
            h4 = h4b + 7.686 ;//btu/lb
            h11 = 1/h1;
            h22 = 1/h2;
            h33 = 1/h3;
            h44 = 1/h4;

            ra = (r*1.8)+32;//in deg f

            kav1 = ra*0.25*(h11+h22+h33+h44);//in fps ---- demand

 //    fvel = (wl*8.33)/(lg*0.07);//in m/s
            fvel_f = (wl*8.33)/(lg*0.07*0.227125);//in fpm conversion factor 0.227125 used

            fh_f = fh/0.3048;
            fs_f = fs*0.00328084;//in ft

            kav2a = 0.06717927* pow(lg,-0.546673);
            kav2b = pow(fvel_f,.25158684);
            kav2c = pow(fh_f,0.6369974);
            kav2d = pow (fs_f,-0.429);
            kav2 = kav2a*kav2b*kav2c*kav2d; //supply ntu

           /* for(lg=1; kav1<=kav2;lg+0.1)
            {
            printf("\n the l/g ratio required for the performance is %f \n",lg);
            printf("\n the KAV/L  for the performance is %f \n",kav2);
            }*/
            /*lg = 1;*/


            //printf("\n the l/g ratio required for the performance is %f \n",lg);
            lg = lg+0.001;
            continue;
               /* if(lg2<lg)
            {
                lg = lg2 + 0.001;

            }
            else{lg=lg+0.2;}*/

            }while(kav1<=kav2 && lg<=3 && lg>=1);
           /* while(lg=1 && lg<=5)
            {
            if(kav2>=kav1)
            break;
            else
            lg = lg+0.2;
            printf("\n the l/g ratio required for the performance is %f \n",lg);
            }*/

               printf("\n the l/g ratio required for the performance is %f \n",lg);
               kav2 = kav2a*kav2b*kav2c*kav2d;
               printf("\n the KAV/L_2  for the performance is %f \n",kav2);
               kav1 = ra*0.25*(h11+h22+h33+h44);
               printf("\n the KAV/L_1  for the performance is %f \n",kav1);

 // fvel = (wl*8.33)/(lg*0.07);//in m/s
               fvel_f = (wl*8.33)/(lg*0.07*0.00508);
               //printf("\n the velocity at the fill is %f m/s or %f fpm \n",fvel,fvel_f);

               l_f = (gpmpc *8.33)/0.227125;//lb/min
               l_s = (gpmpc * 1000)/3600;// kg/sec
               l_f = l_s/0.00755987283;
               printf("\n the mass flow rate of water is %f kg/s or %f lb/min \n",l_s,l_f);

               g_f = l_f/lg;
               g_s = l_s/lg;
               printf("\n the mass flow rate of air is %f kg/s or %f lb/min \n",g_s,g_f);



               printf("\n AT COOLING TOWER INLET\n");

               printf("\n sat vap p inlet is %f Pa",svp_in);
               printf("\n sp humidity inlet is %f kg/kg of dry air",sphu_in);
               spvol_in_s = (spvolro*(tdbi+to)*1000)/(sphum2*(sphupo-svp_in+(sphupo*sphuc*(tdbi-twbi))));//cum/kg
               spvol_in_f = spvol_in_s / 0.0624279606;//cuft/lb
               printf("\n the inlet specific volume is %f cum/kg or %f cu ft/lb ",spvol_in_s,spvol_in_f);//

               row_in_s =(1+ sphu_in)/(spvol_in_s);
               row_in_f = row_in_s/16.0184634;
               printf("\n the inlet density is %f kg/cum or %f lb/cu ft",row_in_s,row_in_f);

               hout_s = spe + (lg*r);
               hout_f = ((hout_s*1.8)/4.1858)+7.686 ;
               printf("\n the inlet enthalpy is %f kj/kg or %f btu/lb",spe_in,hin_f);
               printf("\n the exit enthalpy is %f kj/kg or %f btu/lb",hout_s,hout_f);

               tdbo = hout_s/1.006;
               printf("\n tdbo is %f",tdbo);
           //    tdbo_f = hout_f/0.240;//deg f
               twbo = tdbo;

            too_o = twbo+to;
            svph = svpa/too_o;//
            svpi = svpb*log(too_o);                     //
            svpj = svpd*(twbo-0.01);
            svpk = svpc*(pow(10,svpj));//
            svpl = svpf*(1-(t9/too_o));
            svpm = svpe*(pow(10,svpl));//
            svpn = svph+svpi+svpk+svpm+svpg-0.5 ;
            svp_out = pow(10,svpn);                //sat vap pressure calculated
            printf("\n sat vap p exit is %f Pa   svpn    %f",svp_out,svpn);

            sphua = sphum1/sphum2;//
            sphub = 0;
            sphud = svp_out ;
            sphure = sphud - sphupo ;
            sphuref = abs(sphure);

            sphu_out = (sphua*sphud)/(sphuref);
//          sphuout = (sphua*sphud)/(sphupo - sphud);
            printf("\n sp humidity exit is %f kg/kg of satureated air",sphu_out);


            printf("\n AT COOLING TOWER EXIT\n");

               spvolre = svp_out - sphupo;
            spvolref = abs(spvolre);
               spvol_out_s = (spvolro*(tdbo+to)*1000)/(sphum2*(spvolref));//cum/kg
               spvol_out_f = spvol_out_s / 0.0624279606;//cuft/lb
               printf("\n the exit specific volume is %f cum/kg or %f cu ft/lb ",spvol_out_s,spvol_out_f);

               row_out_s =(1+ sphu_out)/(spvol_out_s);

               row_out_f = row_out_s/16.0184634;
               printf("\n the exit density is %f kg/cum or %f lb/cu ft",row_out_s,row_out_f);

               vel_in_s = (g_s * spvol_in_s)/aia ;//m/s
               vel_in_f = vel_in_s/0.00508;//fpm
               vel_out_s = (g_s*spvol_out_s)/cela ;
               vel_out_f =  vel_out_s/0.00508;

               row_fill_s = (row_in_s + row_out_s)/2 ;
               row_fill_f = (row_in_f + row_out_f)/2 ;
               spvol_fill_s = (spvol_in_s+spvol_out_s)/2;
               spvol_fill_f = (spvol_in_f+spvol_out_f)/2;   //

               vel_fill_s = (g_s*spvol_fill_s)/cela ;
               vel_fill_f = vel_fill_s/0.00508;

               printf("\n the specific volume at fill is %f cum/kg or %f cu ft/lb ",spvol_fill_s,spvol_fill_f);
               printf("\n the density at fill is %f kg/cum or %f lb/cu ft",row_fill_s,row_fill_f);
               printf("\n the velocity of air at inlet is %f m/s or %f fpm ",vel_in_s,vel_in_f);
               printf("\n the velocity of air at fill is %f m/s or %f fpm ",vel_fill_s,vel_fill_f);
               printf("\n the velocity of air at exit is %f m/s or %f fpm ",vel_out_s,vel_out_f);

               cfm_fan_s = g_s * spvol_out_s;//cum/s
               cfm_fan_f = g_f * spvol_out_f;
               printf("\n the volume flow rate of air at fan CFAN is %f cum/s or %f cu ft/min \n",cfm_fan_s,cfm_fan_s);



               pd_il_inh20 = (3*((2*pow(10,-7)*vel_in_f*vel_in_f)+(3*pow(10,-6)*vel_in_f))*(row_in_f/0.07));
        //       pd_il_pa = pd_il_inh20*249;//pascal
               printf("\n pd il %f in ",pd_il_inh20);

               pd_de_inh20 = (5* ((3*pow(10,-7)*vel_out_f*vel_out_f)-(9*pow(10,-6)*vel_out_f) + 0.003 )*(row_out_f/0.07));
        //       pd_de_pa = pd_de_inh20*249;//pascal
               printf("\n pd de %f in ",pd_de_inh20);



               pdf1 = (fh_f/fs_f);
               pdf2 =pow(pdf1,0.963077);
               pdf3 = ((7.44385*pow(10,-6))+((-4.078*pow(10,-6))*pow(fh_f,0.176686)));
               pdf4 = (4.897*(pow(10,-4))*(pow(vel_fill_f,2.3386)));
               pdf5 = (0.0428*wl_f);
               pdf6 = (pow(vel_fill_f,1.33558));

               pd_fill_inh20 = (pdf2*pdf3*(pdf4+(pdf5*pdf6)))*(row_fill_f/0.07);
               printf("\n pd fill %f in ",pd_fill_inh20);
         //      pd_fill_pa = pd_fill_inh20*249;

               tpd_s_inh20 = pd_il_inh20 + pd_de_inh20 + pd_fill_inh20;
               printf("\n hs %f in ",tpd_s_inh20);
               tpd_s_pa = tpd_s_inh20*249;
               pd_v_f = (row_out_f/1202312)*(cfm_fan_f/edisa_f)*(cfm_fan_f/edisa_f);
               printf("\n hv %f in",pd_v_f);
               pd_v_s = pd_v_f * 249 ;
               tpd_f = pd_v_f + tpd_s_inh20;
               tpd_s = tpd_s_pa + pd_v_s;
               printf("\n the overall pressure drop across the cooling tower is %f pa or %f inches of water col",tpd_s,tpd_f);

               bhp = ((cfm_fan_f*(tpd_s_inh20+pd_v_f))/(6356*fann));//hp
               bhp_w = bhp*0.745699872;//kw
               printf("\n the required brake horse power of the fan is %f Bhp or %f kw",bhp,bhp_w);

return (bhp_w);
}

float bestsearch(float delta[vars],float point[vars],float prevbest,int nvars)
{
       float z[vars],minf,ftmp;
       int i;

        minf = prevbest;
for(i = 0; i < nvars ; i++)
{
    z[i] = point[i];
}
    for(i=0 ; i < nvars ; i++)
    {
        z[i] = point[i]+delta[i] ;
        ftmp = bhpeqns(z , nvars);
        if(ftmp < minf)
        {
            minf = ftmp;
        }
        else
        {
            delta[i] = 0.0 - delta[i];
            z[i] = point[i]+delta[i] ;
            ftmp = bhpeqns(z , nvars);
            if(ftmp < minf)
            {
                minf = ftmp;
            }
            else
            {
                z[i] = point[i];
            }
        }
    }
for(i = 0; i < nvars ; i++)
{
     point[i] = z[i];
}

return(minf);
}

float hooke(int nvars,float startpt[vars],float endpt[vars],float rho,float epsilon,int itermax)
{
       float delta[vars];
       float newf, fbefore, steplength, tmp;
       float xbefore[vars], newx[vars];
       int i, j, keep;
       int iters, iadj;
       for (i = 0; i < nvars; i++)
       {
           newx[i] = xbefore[i] = startpt[i];
           delta[i] = fabs(startpt[i] * rho);
           if (delta[i] == 0.0)
           {
               delta[i] = rho;
           }
       }

       iadj = 0;
       steplength = rho;
       iters = 0;
       fbefore = bhpeqns(newx, nvars);
       newf = fbefore;
       while ((iters < itermax) && (steplength > epsilon) && (xbefore[0] >= 0.6096) && (xbefore[0] <= 1.8288) && (xbefore[1] >= 12.0) && (xbefore[1] <= 19.0))
        {
           iters++;
           iadj++;
           printf("\nAfter %5d funevals, f(x) =  %.4le at\n", funevals, fbefore);

           for (j = 0; j < nvars; j++)
           {
               printf("   x[%2d] = %.4le\n", j, xbefore[j]);

           }
              for (i = 0; i < nvars; i++)
           {
               newx[i] = xbefore[i];
           }

           newf = bestsearch(delta, newx, fbefore, nvars);
           keep = 1;
              while ((newf < fbefore) && (keep == 1))
             {
               iadj = 0;
               for (i = 0; i < nvars; i++)
                 {
                   if(newx[i] <= xbefore[i])
                    {
                     delta[i] = 0.0 - fabs(delta[i]);
                    }
                    else
                    {
                       delta[i] = fabs(delta[i]);
                    }
                   tmp = xbefore[i];
                   xbefore[i] = newx[i];
                   newx[i] = newx[i] + newx[i] - tmp;
               }

               fbefore = newf;
               newf = bestsearch(delta, newx, fbefore, nvars);

               if (newf >= fbefore)
                {
                   break;
                 }
               keep = 0;
               for (i = 0; i < nvars; i++)
                {
                   keep = 1;
                   if (fabs(newx[i] - xbefore[i]) > (0.5 * fabs(delta[i])))
                    {
                      break;
                     }
                   else
                     {
                       keep = 0;
                     }
               }


             if ((steplength >= epsilon) && (newf >= fbefore))
             {
               steplength = steplength * rho;
               for (i = 0; i < nvars; i++)
               {
                   delta[i] *= rho;
               }
             }
       }
    for (i = 0;(i < nvars) ; i++)
    {
          endpt[i] = xbefore[i];
    }

}
 getch();
 return(iters);
}

this program is running but not giving the desired results and i didnt use the calloc and exit codes as my compiler is throwing an error there....the result for x[0] and x[1] is not in the desired range of the values

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 10 · 2013-04-11T15:25:29+00:00

If you don't mind me asking, which compiler are you using? If it is Turbo C++ (or any other 16-bit compiler), there may be issues with the size of the int and float variables being used. If I recall correctly, Turbo C uses 16-bit int values and 32-bit floats.

sai.aleenus 0 Newbie Poster · Answer 11 · 2013-04-11T18:16:09+00:00

sai.aleenus 0 Newbie Poster

11 Years Ago

i am using dosbox...

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 12 · 2013-04-11T20:57:22+00:00

Ah, DosBox isn't a compiler, it is a utility that allows you to run 16-bit programs (e.g., Turbo C++) in a modern Windows environment. I will assume that it is the Turbo compiler, then, or something from the same time period (c. 1990).

This could present a problem, as it is possible that you are getting an overflow (or an underflow) somewhere, and it isn't getting detected. I don't know if that's the case, but it may be, which case I'm not sure what you can do about it other than getting a newer compiler.

BTW, is this supposed to be in C, or C++? You use C++ style comments in several places, which most C compilers won't accept unless they support the newest language standard; yet you aren't using any C++ functions or operators. Just wondering.

Addressing the code isues, I again want to stress dividing the bhpeqns() function up into smaller functions. As it is, it sprawls over 370 LOC, and contains well over 100 discrete variables. This is simply unmanageable. If you want any chance of getting this to work, start breaking the sections of code up into functions, like so:

float saturation_vapor_pressure(float vapor[], float twbi, float to, float t9)
{
    float too, saturation[7];

    too = twbi + to;
    saturation[0] = vapor[0] / too;
    saturation[1] = vapor[1] * log(too);
    saturation[2] = vapor[2] * (twbi - 0.01);
    saturation[3] = vapor[3] * pow(10, saturation[3]);
    saturation[4] = vapor[4] * (1 - (t9 / too));
    saturation[5] = vapor[5] * pow(10, saturation[4]);//
    saturation[6] = vapor[6] + saturation[0] + saturation[1] + saturation[3] + saturation[5];

    return pow(10, saturation[6]);                //sat vap pressure calculated
}

I specifically chose this function in part because it is one you repeat later in the code, with separate variables; both of the instances can now be replaced by a single function call each.

I'm guessing with the names of the variables; you may have better names for them yourself, and if you do I would recommend using them. If it would help any, gather related values up into arrays or structures wherever possible.

One more question: are the initialized values starting at line 25 variables, or constants? If they aren't changed anywhere in the code, you should declare them as constants, or else use #define to name them.

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 13 · 2013-04-11T21:34:20+00:00

gather related values up into arrays or structures wherever possible.

I forgot to give an example of what I meant by this. What I meant was, that instead of something like this:

float svpa = -2948.997118, svpb = -2.1836674, svpc = -0.000150474,svpd = -0.0303738468;
float svpe = 0.00042873, svpf = 4.76955,svpg = 25.83220018;

use this:

float vapor[] = {-2948.997118, -2.1836674, -0.000150474, -0.0303738468, 0.00042873, 4.76955, 5.83220018};

or this:

struct VAPOR 
{
    float a, b, c, d, e, f, g;
} vapor = {-2948.997118, -2.1836674, -0.000150474, -0.0303738468, 0.00042873, 4.76955, 5.83220018};

It helps organize the values a little better.

sai.aleenus 0 Newbie Poster · Answer 14 · 2013-04-12T03:24:16+00:00

i will make the function use as suggested by you..
the initialised values are not constant and ther were only one set of data...i am supposed to vary those values as well and note the variations..

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 15 · 2013-04-12T11:49:46+00:00

When you say you were expected to vary the results, do you mean you would set them in the program differently and recompile and run it, or that you would need to be able to enter in the data from somewhere (presumably a data file - entering all that in by hand would be certain to result in frequent errors) and pass the data to the functions? I ask this because as things are, you are basically treating the initialized values as constants - they never get reassigned at any point - which is why I thought they might be constants in the program.

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 16 · 2013-04-12T16:21:42+00:00

I've been working on this most of the morning, and have at least re-organized some aspects of it into something that should be reasonable; unfortunately, the results I get are anything but. One of the figures is overflowing right at the start, and this is throwing off the whole program. Perhaps you can use what I've done to make the program fit the data better.

hooke-jeeves-test.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdbool.h>
#include "hooke-jeeves.h"


#define VARS (250)
#define RHO (0.85)
#define EPSILON (0.0001)
#define ITERMAX (1024)
#define NVARS (2)


bool initialize_plant_state(FILE* data_src);

double bhpeqns(double x[], int n, const int vars);

double saturation_vapor_pressure(double twbi);
double sp_humidity(double d, double w, double svp);
double enthalpy_metric(double d, double humid);
double enthalpy_btu(double enthalpy);


// global 'constants' to be set in initialize_plant_state() -
// they vary with the initial data set, but are not altered during the
// program run

double thw, tcw, twbi_init, tdbi, fand, hubd, fann;
double aih, cell, celw, gpm;
int ais, celn;
double to, t9;
double vapor[7];
struct
{
    double hum1, hum2, po, c;
} humidity;

double spec2, spec3, spelo, spvolro;
int wlmin, wlmax, vmin, vmax;


// global variables - these are set in the functions
double app, r;
double aia, edisa, cela, gpmpc, wl;


int main()
{
    double startpt[VARS], endpt[VARS];
    int i, jj;

    char filename[FILENAME_MAX];
    FILE* data_src;
    bool success;

    printf("Enter the name of the data file: ");
    scanf("%4095s", filename);
    data_src = fopen(filename, "r");

    if (data_src == NULL)
    {
        printf("Unable to open file '%s', exiting.\n", filename);
        exit(1);
    }

    success = initialize_plant_state(data_src);
    fclose(data_src);

    if (!success)
    {
        printf("Could not read all the require data, exiting.\n");
        exit(1);
    }


    /* starting guess for rosenbrock test function */
    startpt[0] = 0.6892;
    startpt[1] = 12.0f;
    jj = hooke(bhpeqns, VARS, NVARS, startpt, endpt, RHO, EPSILON, ITERMAX);
    printf("\n\n\nHOOKE USED %d ITERATIONS, AND RETURNED\n", jj);
    for (i = 0; i < NVARS; i++)
    {
        printf("x[%3d] = %15.7e \n", i, endpt[i]);
    }

    return 0;
}



bool initialize_plant_state(FILE* data_src)
{
    if (7 > fscanf(data_src, "%lf %lf %lf %lf %lf %lf %lf", &thw, &tcw, &twbi_init, &tdbi, &fand, &hubd, &fann))
    {
        printf("line 0: could not read one or more of the following: %lf %lf %lf %lf %lf %lf %lf\n", thw, tcw, twbi_init, tdbi, fand, hubd, fann);
        return false;
    }

    if (4 > fscanf(data_src, "%lf %lf %lf %lf\n", &aih, &cell, &celw, &gpm))
    {
        printf("line 1: could not read one or more of the following: %lf %lf %lf %lf\n", aih, cell, celw, gpm);
        return false;
    }
    if (2 > fscanf(data_src, "%d %d\n", &ais, &celn))
    {
        printf("line 2: could not read one or more of the following: %d %d\n", ais, celn);
        return false;
    }
    if (2 > fscanf(data_src, "%lf %lf\n", &to, &t9))
    {
        printf("line 3: could not read one or more of the following: %d %d\n", to, t9);
        return false;
    }
    if (7 > fscanf(data_src, "%lf %lf %lf %lf %lf %lf %lf\n", &vapor[0], &vapor[1], &vapor[2], &vapor[3], &vapor[4], &vapor[5], &vapor[6]))
    {
        printf("line 4: could not read one or more of the following: %lf %lf %lf %lf %lf %lf %lf\n",
                    vapor[0], vapor[1], vapor[2], vapor[3], vapor[4], vapor[5], vapor[6]);
        return false;
    }
    if (4 > fscanf(data_src, "%lf %lf %lf %lf\n", &humidity.hum1, &humidity.hum2, &humidity.po, &humidity.c))
    {
        printf("line 5: could not read one or more of the following: %lf %lf %lf %lf\n", humidity.hum1, humidity.hum2, humidity.po, humidity.c);
        return false;
    }
    if (4 > fscanf(data_src, "%lf %lf %lf %lf\n", &spec2, &spec3, &spelo, &spvolro))
    {
        printf("line 6: could not read one or more of the following: %lf %lf %lf %lf\n", spec2, spec3, spelo, spvolro);
        return false;
    }
    if (4 > fscanf(data_src, "%d %d %d %d\n", &wlmin, &wlmax, &vmin, &vmax))
    {
        printf("line 7: could not read one or more of the following: %d %d %d %d\n", wlmin, wlmax, vmin, vmax);
        return false;
    }

    return true;
}


double saturation_vapor_pressure(double twbi)
{
    double too, saturation[7];

    too = twbi + to;
    saturation[0] = vapor[0] / too;
    saturation[1] = vapor[1] * log(too);
    saturation[2] = vapor[2] * (twbi - 0.01);
    saturation[3] = vapor[3] * pow(10, saturation[3]);
    saturation[4] = vapor[4] * (1 - (t9 / too));
    saturation[5] = vapor[5] * pow(10, saturation[4]);//
    saturation[6] = vapor[6] + saturation[0] + saturation[1] + saturation[3] + saturation[5];

    return pow(10, saturation[6]);                //sat vap pressure calculated
}



double sp_humidity(double d, double w, double svp)
{
    double sphua, sphub, sphud;

    sphua = humidity.hum1 / humidity.hum2;
    sphub = humidity.po * humidity.c * (d - w);
    sphud = svp - sphub;
    return (sphua * sphud) / (humidity.po - sphud); // sp humidity is calculated
}


double enthalpy_metric(double d, double humid)
{
    return (spec2 * d) + (humid * (spelo + (spec3 * d)));
}

double enthalpy_btu(double enthalpy)
{
    return ((enthalpy * 1.8) / 4.1858) + 7.686;
}


double bhpeqns(double x[], int n, const int vars)
{
    int i;

    double fh_f,fs_f,fh,fs;
    double svp[5];
    double hin_f[5];
    double twbi[5], t[5];
    double sphu[5],spe[5];
    double h[4], h_inv[4], h_inv_sum = 0.0;
    double speat[4], ra, ha[4], hb[4];
    double kav1,kav2,lg;
    double kav2a,kav2b,kav2c,kav2d,edisa_f;
    double fvel_f,l_f,l_s,g_f,g_s,spvol_in_s,spvol_in_f,row_in_s,row_in_f;
    double twbo,tdbo,svp_out,sphu_out,spvol_out_s,spvol_out_f,hout_s,hout_f,row_out_s,row_out_f;
    double vel_in_s,vel_in_f,vel_out_s,vel_out_f,vel_fill_s,vel_fill_f,row_fill_s,row_fill_f,spvol_fill_s,spvol_fill_f;
    double wl_f;
    double pdf[6];
    double pd_il_inh20,cfm_fan_s,cfm_fan_f,pd_de_inh20, pd_fill_inh20;
    double tpd_s_inh20,tpd_s_pa,pd_v_f,pd_v_s,bhp,bhp_w,tpd_f,tpd_s;
    double spvolref,spvolre;
    double edisa_s;


    printf("\n %d",n);

    t[0] = tcw;
    t[1] = tcw + (0.1 * r);
    t[2] = tcw + (0.4 * r);
    t[3] = thw - (0.4 * r);
    t[4] = thw - (0.1 * r);

    twbi[0] = twbi_init;
    twbi[1] = twbi[0] + (0.1 * r);
    twbi[2] = twbi[0] + (0.4 * r);
    twbi[3] = twbi[0] + (0.6 * r);
    twbi[4] = twbi[0] + (0.9 * r);

    edisa_s = pow(3.2808, 2);

    r = thw - tcw;
    app = tcw - twbi[0] ;
    cela = cell * celw;
    aia = ais * celw * aih;
    edisa = 0.785 * (pow(fand, 2) - pow(hubd, 2));
    edisa_f = edisa * edisa_s;
    printf("\n edisa %lf",edisa_f);
    gpmpc = gpm / celn ;
    wl = gpmpc / cela;
    wl_f = (gpmpc * 4.402867) / (cela * edisa_s);// fpm
    fh = x[0];
    fs = x[1];


    for (i = 0; i < 5; i++)
    {

        svp[i] = pow(10, saturation_vapor_pressure(twbi[i]));                //sat vap pressure calculated
        printf("\n sat vap pressure is %lf Pa", svp[i]);

        sphu[i] = sp_humidity(tdbi, twbi[i], svp[i]);
        printf("\n sp humidity is %lf kg/kg of dry air",sphu[i]);

        spe[i] = enthalpy_metric(tdbi, sphu[i]);//sp enthalpy is calculated
        printf("\n enthalpy is(kj/kg) %lf",spe[i]);

        hin_f[i] = enthalpy_btu(spe[i]);
        printf("\n enthalpy is(btu/lb) %lf", hin_f[i]);
    }


    lg = 1.0f;

    do
    {
        speat[0] = spe[4] + (0.1 * lg * r);
        speat[1] = spe[4] + (0.4 * lg * r);
        speat[2] = spe[4] + (0.6 * lg * r);
        speat[3] = spe[4] + (0.9 * lg * r);

        for (i = 0; i < 4; i++)
        {
            ha[i] = spe[i + 1] - speat[i];
            hb[i] = ((ha[i] * 1.8) / 4.1858);
            h[i] = hb[i] + 7.686 ;
            h_inv[i] = 1 / h[i];
            h_inv_sum += h_inv[i];
        }

        ra = (r*1.8)+32;//in deg f

        kav1 = ra * 0.25 * (h_inv_sum);//in fps ---- demand

        fvel_f = (wl * 8.33) / (lg * 0.07 * 0.227125);  //in fpm conversion factor 0.227125 used

        fh_f = fh / 0.3048;
        fs_f = fs * 0.00328084;//in ft

        kav2a = 0.06717927 * pow(lg, -0.546673);
        kav2b = pow(fvel_f, 0.25158684);
        kav2c = pow(fh_f, 0.6369974);
        kav2d = pow (fs_f, -0.429);
        kav2 = kav2a * kav2b * kav2c * kav2d; //supply ntu

        lg += 0.001;
    }
    while(kav1<=kav2 && lg<=3 && lg>=1);


    printf("\n the l/g ratio required for the performance is %lf \n",lg);
    kav2 = kav2a*kav2b*kav2c*kav2d;
    printf("\n the KAV/L_2  for the performance is %lf \n",kav2);
    kav1 = ra * 0.25 * h_inv_sum;
    printf("\n the KAV/L_1  for the performance is %lf \n",kav1);


    fvel_f = (wl * 8.33) / (lg * 0.07 * 0.00508);


    l_f = (gpmpc * 8.33) / 0.227125;//lb/min
    l_s = (gpmpc * 1000) / 3600;// kg/sec
    l_f = l_s / 0.00755987283;
    printf("\n the mass flow rate of water is %lf kg/s or %lf lb/min \n",l_s,l_f);

    g_f = l_f / lg;
    g_s = l_s / lg;
    printf("\n the mass flow rate of air is %lf kg/s or %lf lb/min \n",g_s,g_f);



    printf("\n AT COOLING TOWER INLET\n");

    printf("\n sat vap p inlet is %lf Pa",svp[0]);
    printf("\n sp humidity inlet is %lf kg/kg of dry air", sphu[0]);
    spvol_in_s = (spvolro * (tdbi + to) * 1000) / (humidity.hum2 * (humidity.po - svp[0] + (humidity.po * humidity.c * (tdbi - twbi[0]))));//cum/kg
    spvol_in_f = spvol_in_s / 0.0624279606;//cuft/lb
    printf("\n the inlet specific volume is %lf cum/kg or %lf cu ft/lb ",spvol_in_s,spvol_in_f);//

    row_in_s =(1+ sphu[0])/(spvol_in_s);
    row_in_f = row_in_s/16.0184634;
    printf("\n the inlet density is %lf kg/cum or %lf lb/cu ft",row_in_s,row_in_f);

    hout_s = spe[5] + (lg*r);
    hout_f = ((hout_s * 1.8) / 4.1858) + 7.686 ;
    printf("\n the inlet enthalpy is %lf kj/kg or %lf btu/lb", spe[0], hin_f[0]);
    printf("\n the exit enthalpy is %lf kj/kg or %lf btu/lb", hout_s, hout_f);

    tdbo = hout_s / 1.006;
    printf("\n tdbo is %lf",tdbo);
    //    tdbo_f = hout_f/0.240;//deg f
    twbo = tdbo;

    svp_out = pow(10, saturation_vapor_pressure(twbo));                //sat vap pressure calculated
    printf("\n sat vap p exit is %lf Pa   svpn    %lf",svp_out, svp[0]);

    sphu_out = sp_humidity(0, 0, svp_out);
//          sphuout = (sphua*sphud)/(sphupo - sphud);
    printf("\n sp humidity exit is %lf kg/kg of saturated air", sphu_out);


    printf("\n AT COOLING TOWER EXIT\n");

    spvolre = svp_out - humidity.po;
    spvolref = abs(spvolre);
    spvol_out_s = (spvolro*(tdbo+to)*1000)/(humidity.hum2*(spvolref));//cum/kg
    spvol_out_f = spvol_out_s / 0.0624279606;//cuft/lb
    printf("\n the exit specific volume is %lf cum/kg or %lf cu ft/lb ",spvol_out_s,spvol_out_f);

    row_out_s =(1+ sphu_out)/(spvol_out_s);

    row_out_f = row_out_s/16.0184634;
    printf("\n the exit density is %lf kg/cum or %lf lb/cu ft",row_out_s,row_out_f);

    vel_in_s = (g_s * spvol_in_s)/aia ;//m/s
    vel_in_f = vel_in_s/0.00508;//fpm
    vel_out_s = (g_s*spvol_out_s)/cela ;
    vel_out_f =  vel_out_s/0.00508;

    row_fill_s = (row_in_s + row_out_s)/2 ;
    row_fill_f = (row_in_f + row_out_f)/2 ;
    spvol_fill_s = (spvol_in_s+spvol_out_s)/2;
    spvol_fill_f = (spvol_in_f+spvol_out_f)/2;   //

    vel_fill_s = (g_s*spvol_fill_s)/cela ;
    vel_fill_f = vel_fill_s/0.00508;

    printf("\n the specific volume at fill is %lf cum/kg or %lf cu ft/lb ",spvol_fill_s,spvol_fill_f);
    printf("\n the density at fill is %lf kg/cum or %lf lb/cu ft",row_fill_s,row_fill_f);
    printf("\n the velocity of air at inlet is %lf m/s or %lf fpm ",vel_in_s,vel_in_f);
    printf("\n the velocity of air at fill is %lf m/s or %lf fpm ",vel_fill_s,vel_fill_f);
    printf("\n the velocity of air at exit is %lf m/s or %lf fpm ",vel_out_s,vel_out_f);

    cfm_fan_s = g_s * spvol_out_s;//cum/s
    cfm_fan_f = g_f * spvol_out_f;
    printf("\n the volume flow rate of air at fan CFAN is %lf cum/s or %lf cu ft/min \n",cfm_fan_s,cfm_fan_s);



    pd_il_inh20 = (3*((2*pow(10,-7)*vel_in_f*vel_in_f)+(3*pow(10,-6)*vel_in_f))*(row_in_f/0.07));
    //       pd_il_pa = pd_il_inh20*249;//pascal
    printf("\n pd il %lf in ",pd_il_inh20);

    pd_de_inh20 = (5* ((3*pow(10,-7)*vel_out_f*vel_out_f)-(9*pow(10,-6)*vel_out_f) + 0.003 )*(row_out_f/0.07));
    //       pd_de_pa = pd_de_inh20*249;//pascal
    printf("\n pd de %lf in ",pd_de_inh20);



    pdf[0] = (fh_f/fs_f);
    pdf[1] = pow(pdf[0], 0.963077);
    pdf[2] = ((7.44385*pow(10,-6))+((-4.078*pow(10,-6))*pow(fh_f,0.176686)));
    pdf[3] = (4.897*(pow(10,-4))*(pow(vel_fill_f,2.3386)));
    pdf[4] = (0.0428*wl_f);
    pdf[5] = (pow(vel_fill_f,1.33558));

    pd_fill_inh20 = (pdf[1] * pdf[2] * (pdf[3] + (pdf[4] * pdf[5]))) * (row_fill_f / 0.07);
    printf("\n pd fill %lf in ",pd_fill_inh20);
    //      pd_fill_pa = pd_fill_inh20*249;

    tpd_s_inh20 = pd_il_inh20 + pd_de_inh20 + pd_fill_inh20;
    printf("\n hs %lf in ",tpd_s_inh20);
    tpd_s_pa = tpd_s_inh20*249;
    pd_v_f = (row_out_f/1202312)*(cfm_fan_f/edisa_f)*(cfm_fan_f/edisa_f);
    printf("\n hv %lf in",pd_v_f);
    pd_v_s = pd_v_f * 249 ;
    tpd_f = pd_v_f + tpd_s_inh20;
    tpd_s = tpd_s_pa + pd_v_s;
    printf("\n the overall pressure drop across the cooling tower is %lf pa or %lf inches of water col",tpd_s,tpd_f);

    bhp = ((cfm_fan_f*(tpd_s_inh20+pd_v_f))/(6356*fann));//hp
    bhp_w = bhp*0.745699872;//kw
    printf("\n the required brake horse power of the fan is %lf Bhp or %lf kw",bhp,bhp_w);

    return (bhp_w);
}

The data file is as follows:

43.0 40.0 37.0 40.0 7.0 1.0 0.75
3.0 8.0 8.0 1000.0
3 4
273.15 273.16
-2948.997118 -2.1836674 -0.000150474 -0.0303738468 0.00042873 4.76955 25.83220018
18.01534 28.9645 101325 0.000666
1.00568 1.84598 2500.84 8.31432
4 11 300 750

I'd like to hear back from you about this, to see if you can shed any light on the unusual results (or at least give me some idea of what the variables mean); in the meanwhile, I need to work on other projects right now, but I'll get back to you sometime today.

sai.aleenus 0 Newbie Poster · Answer 17 · 2013-04-12T18:29:41+00:00

my basic program is c and not cpp..i have followed your suggestion and grouped the data properly into structures and have triedthe code again..i am not encountering any problem with the bhpeqns function but am encountering the problem at hooke..i am attaching the code i have done till far...
i am also giving the step wise priper description of my project so as to throw light on what exactly i am doing...
aim: to optimise the bhpeqns function of my program(hkjv15) using hooke and jeeves optimisation.

imp : the bhp function has 2 arguments fh and fs ,which have ranges : 0.6096<=fh<=1.8288 and 12<=fs<=19 and the final values of fh and fs obtained for the optimum bhp should lie between their given ranges.

the aim is to find the values of fh and fs where the function bhp value is the minimum and also the values of fh and fs lie in their corresponding ranges.

*****logic of hooke and jeeve:***

(note initially delta is a const of 0.8 for fh and 1 for fs)

1.initially a value of fh and fs are inputted within their ranges and initial bhp (bhpi) is calculated by bhpeqns(fh,fs) and the ranges of fh and fs are to be implemented.

2.x direction move is implemented:

fh1 = fh + (delta * 1) or fh1 = fh - (delta * 1) while fs1 = fs

we determine the + or – sign of fh by the one that yields the bhp<bhpi

bhp1 = bhpeqns(fh1,fs1) note that delta is to be halved till bhp1<bhpi

3.y direction move is implemented:

fh2 = fh1 while fs2 = fs1 +/-(delta *1)

we determine the + or – sign of fs by the one that yields the bhp2<bhp1

bhp2 = bhpeqns(fh2,fs2)

4.pattern search

fh pattern directon fhsp = fh2-fh1

fs pattern direction fssp = fs2-fs1

rho is a const = lam

fh3 = fh2 + (fhsp * lam)

fs3 = fs2 + (fssp * lam)

bhp3 = bhpeqns(fh3,fs3)

5.if bhp3<bhp2 steps 2,3,4 are to be repeated till the lowest bhp is obtained with fh = fh3 and fs = fs2

else fh = fh2 and fs = fs2 and the steps 2,3,4 are to be repeated.

i am running my code but it is running an infinite loop though i am curtailing the number of iterations.
please help me with the code in the function of hooke

#include <stdio.h>
#include <stdlib.h>
#include<conio.h>
#include<math.h>
#include<float.h>
#include<stdarg.h>
//#include<graphics.h>
#include<iostream.h>;

#define to 273.15
#define t9 273.16
#define spvolro 8.31432

float sat_vap(float twbi);
float sp_hum(float svp,float tdbi,float twbi);
float sat_vape(float twbo);
float sp_enthalpy(float tdbi,float sphum);
float sp_vol(float tdbi,float twbi,float svp);
float row(float sphum,float spvol);

float bhpeqns(float fh,float fs);

float fsdet(float fs);

float fhdet(float fh)
float hooke(int max,float fh,float fs,float epsilon,float fhsl);



float thw=43.0f,tcw=40.0f,twbi=37.0f,tdbi=40.0f,app,r,fand=7.0f,hubd=1.0f,fann=0.75f;
float aia,aih=3.0f,edisa,edisa_f,cell=8.0f,celw=8.0f,cela,twbo,tdbo,gpm=1000.0f,gpmpc,wl,spe_in,wl_f;
int ais = 3,celn = 4;
float svpp[] = { -2948.997118,  -2.1836674, -0.000150474, -0.0303738468,0.00042873,  4.76955, 25.83220018,-0.5};
float sphup[] = {18.01534, 28.9645,  101325,  0.000666};    //m1,m2,po,c
float spep[] = {  1.00568, 1.84598, 2500.84};       //c2,c3,lo
float  spec2 = 1.00568, spec3 = 1.84598, spelo = 2500.84,spvolro = 8.31432;float row_fill_f;
int i,j,wlmin=4 , wlmax = 11, vmin = 300, vmax = 750;
int k;
// float wl_f,pdf1,pdf2,pdf3,pdf4,pdf5,pdf6;



#define fhsl_ini 0.85
#define epsimin 0.0001
#define imax 100






void main()
{
  float fh,fs,fhsl,epsilon;
  int ima,k;
  printf("\n enter a value for fill height(note in b/w 0.6096 and 1.8288) : \t");
  scanf("%f",&fh);
  printf("\n enter a value for flute spacing(note in b/w 12 and 19) : \t");
  scanf("%f",&fs);
  fh = fh;
  fs = fs;
  fhsl = fhsl_ini;
  ima = imax;
  epsilon = epsimin;
  k = hooke(ima,fh,fs,epsilon,fhsl);
  printf("\n the iterations done are %d",k);
getch();


}


float bhpeqns(float fh, float fs)
{

    printf("i got called");

    float fh_f,fs_f;
    float svp_in,svph,svpi,svpj,svpk, svpl, svpm, svpn,too;//sat vap pressure
    float sphu_in,sphua,sphub,sphud,hin_f;//sp humidity
    float too1,too2,too3,too4,twbi1,twbi2,twbi3,twbi4,svp1,svp2,svp3,svp4,t1,t2,t3,t4,tdbi1,tdbi2,tdbi3,tdbi4;
    float sphu1,sphu2,sphu3,sphu4,spe1,spe2,spe3,spe4,spe;
    float h1,h2,h3,h4,h11,h22,h33,h44;
    float speat1,speat2,speat3,speat4,ra,h1a,h1b,h2a,h2b,h3a,h3b,h4a,h4b;
    float kav1,kav2,lg;
    float kav2a,kav2b,kav2c,kav2d,edisa_f;
    float fvel,fvel_f,l_f,l_s,g_f,g_s,spvol_in_s,spvol_in_f,row_in_s,row_in_f;
    float too_o,twbo,tdbo,svp_out,sphu_out,spvol_out_s,spvol_out_f,hout_s,hout_f,tdbo_f,row_out_s,row_out_f;
    float vel_in_s,vel_in_f,vel_out_s,vel_out_f,vel_fill_s,vel_fill_f,row_fill_s,row_fill_f,spvol_fill_s,spvol_fill_f;
    float wl_f,pdf1,pdf2,pdf3,pdf4,pdf5,pdf6,pd_il_inh20,pd_il_pa,cfm_fan_s,cfm_fan_f,pd_de_inh20,pd_de_pa,pd_fill_inh20,pd_fill_pa;
    float tpd_s_inh20,tpd_s_pa,pd_v_f,pd_v_s,bhp,bhp_w,tpd_f,tpd_s;
    float sphure,sphuref,spvolref,spvolre,predp;
    float pdvf2,pdvf3,bhhp1,bhhp2;



            r = thw - tcw ;
            app = tcw - twbi ;
            cela = cell*celw ;
            aia = ais*celw*aih;
            edisa = (0.785*((fand*fand)-(hubd*hubd)));
            edisa_f = edisa*3.2808*3.2808;
            printf("\n edisa %f",edisa_f);
            gpmpc = gpm/celn ;
            wl = gpmpc/cela;
            wl_f = (gpmpc*4.402867)/(cela*3.2808*3.2808);// fpm


            twbi = twbi;    
            svp_in = sat_vap(twbi);                //sat vap pressure calculated
            printf("\n sat vap p inlet is %f Pa",svp_in);

            tdbi = tdbi;
            twbi = twbi;    
            sphu_in = sp_hum(svp_in, tdbi, twbi); // sp humidity is calculated
            printf("\n sp humidity inlet is %f kg/kg of dry air",sphu_in);

            spe_in = sp_enthalpy(tdbi,sphu_in);//sp enthalpy is calculated
            printf("\n inlet enthalpy is(kj/kg) %f",spe_in);
            spe = spe_in;

            hin_f = ((spe_in*1.8)/4.1858)+7.686 ;
            printf("\n inlet enthalpy is(btu/lb) %f",hin_f);

            lg = 1.0f;
            do
            {

            t1 = tcw + (0.1*r);
            t2 = tcw + (0.4*r);
            t3 = thw - (0.4*r);
            t4 = thw - (0.1*r);

            twbi1 = twbi + (0.1*r);
            twbi2 = twbi + (0.4*r);
            twbi3 = twbi + (0.6*r);
            twbi4 = twbi + (0.9*r);

            tdbi1 = twbi + (0.1*r);
            tdbi2 = twbi + (0.4*r);
            tdbi3 = twbi + (0.6*r);
            tdbi4 = twbi + (0.9*r); 

            svp1 = sat_vap(twbi1);                //sat vap pressure calculated
         //   printf("\n sat vap p1 is %f",svp1);


            sphu1 = sp_hum(svp1, tdbi1, twbi1);  // sp humidity is calculated
          //  printf("\n sp humidity1 is %f",sphu1);


            spe1 = sp_enthalpy(tdbi1,sphu1);//sp enthalpy is calculated
          //  printf("\n enthalpy1 is %f",spe1);

            svp2 = sat_vap(twbi2);                //sat vap pressure calculated
         //   printf("\n sat vap p2 is %f",svp2);


            sphu2 = sp_hum(svp2, tdbi2, twbi2);  // sp humidity is calculated
          //  printf("\n sp humidity3 is %f",sphu2);


            spe2 = sp_enthalpy(tdbi2,sphu2);//sp enthalpy is calculated
          //  printf("\n enthalpy2 is %f",spe2);



            svp3 = sat_vap(twbi3);                //sat vap pressure calculated
         //   printf("\n sat vap p3 is %f",svp3);


            sphu3 = sp_hum(svp3, tdbi3, twbi3);  // sp humidity is calculated
          //  printf("\n sp humidity3 is %f",sphu3);


            spe3 = sp_enthalpy(tdbi3,sphu3);//sp enthalpy is calculated
          //  printf("\n enthalpy3 is %f",spe1);


            svp4 = sat_vap(twbi4);                //sat vap pressure calculated
         //   printf("\n sat vap p4 is %f",svp4);


            sphu4 = sp_hum(svp4, tdbi4, twbi4);  // sp humidity is calculated
          //  printf("\n sp humidity4 is %f",sphu4);


            spe4 = sp_enthalpy(tdbi4,sphu4);//sp enthalpy is calculated
          //  printf("\n enthalpy4 is %f",spe4);

            speat1 = spe + (0.1*lg*r);
            speat2 = spe + (0.4*lg*r);
            speat3 = spe + (0.6*lg*r);
            speat4 = spe + (0.9*lg*r);

            h1a = spe1-speat1;
            h1b = ((h1a*1.8)/4.1858);
            h1 = h1b + 7.686 ;
            h2a = spe2-speat2;
            h2b = ((h2a*1.8)/4.1858);
            h2 = h2b + 7.686 ;
            h3a = spe3-speat3;
            h3b = ((h3a*1.8)/4.1858);
            h3 = h3b + 7.686 ;
            h4a = spe4-speat4;
            h4b = ((h4a*1.8)/4.1858);
            h4 = h4b + 7.686 ;//btu/lb
            h11 = 1/h1;
            h22 = 1/h2;
            h33 = 1/h3;
            h44 = 1/h4;

            ra = (r*1.8)+32;//in deg f

            kav1 = ra*0.25*(h11+h22+h33+h44);//in fps ---- demand

 //    fvel = (wl*8.33)/(lg*0.07);//in m/s
            fvel_f = (wl*8.33)/(lg*0.07*0.227125);//in fpm conversion factor 0.227125 used

            fh_f = fh/0.3048;
            fs_f = fs*0.00328084;//in ft

            kav2a = 0.06717927* pow(lg,-0.546673);
            kav2b = pow(fvel_f,.25158684);
            kav2c = pow(fh_f,0.6369974);
            kav2d = pow (fs_f,-0.429);
            kav2 = kav2a*kav2b*kav2c*kav2d; //supply ntu

           /* for(lg=1; kav1<=kav2;lg+0.1)
            {
            printf("\n the l/g ratio required for the performance is %f \n",lg);
            printf("\n the KAV/L  for the performance is %f \n",kav2);
            }*/
            /*lg = 1;*/


            //printf("\n the l/g ratio required for the performance is %f \n",lg);
            lg = lg+0.001;
            continue;
               /* if(lg2<lg)
            {
                lg = lg2 + 0.001;

            }
            else{lg=lg+0.2;}*/

            }while(kav1<=kav2 && lg<=3 && lg>=1);
           /* while(lg=1 && lg<=5)
            {
            if(kav2>=kav1)
            break;
            else
            lg = lg+0.2;
            printf("\n the l/g ratio required for the performance is %f \n",lg);
            }*/

               printf("\n the l/g ratio required for the performance is %f \n",lg);
               kav2 = kav2a*kav2b*kav2c*kav2d;
               printf("\n the KAV/L_2  for the performance is %f \n",kav2);
               kav1 = ra*0.25*(h11+h22+h33+h44);
               printf("\n the KAV/L_1  for the performance is %f \n",kav1);

 // fvel = (wl*8.33)/(lg*0.07);//in m/s
               fvel_f = (wl*8.33)/(lg*0.07*0.00508);
               //printf("\n the velocity at the fill is %f m/s or %f fpm \n",fvel,fvel_f);

               l_f = (gpmpc *8.33)/0.227125;//lb/min
               l_s = (gpmpc * 1000)/3600;// kg/sec
               l_f = l_s/0.00755987283;
               printf("\n the mass flow rate of water is %f kg/s or %f lb/min \n",l_s,l_f);

               g_f = l_f/lg;
               g_s = l_s/lg;
               printf("\n the mass flow rate of air is %f kg/s or %f lb/min \n",g_s,g_f);



               printf("\n AT COOLING TOWER INLET\n");

               printf("\n sat vap p inlet is %f Pa",svp_in);
               printf("\n sp humidity inlet is %f kg/kg of dry air",sphu_in);
               spvol_in_s = sp_vol(tdbi,twbi, svp_in);;//cum/kg
               spvol_in_f = spvol_in_s / 0.0624279606;//cuft/lb
               printf("\n the inlet specific volume is %f cum/kg or %f cu ft/lb ",spvol_in_s,spvol_in_f);//

               row_in_s = row(sphu_in,spvol_in_s);;
               row_in_f = row_in_s/16.0184634;
               printf("\n the inlet density is %f kg/cum or %f lb/cu ft",row_in_s,row_in_f);

               hout_s = spe + (lg*r);
               hout_f = ((hout_s*1.8)/4.1858)+7.686 ;
               printf("\n the inlet enthalpy is %f kj/kg or %f btu/lb",spe_in,hin_f);
               printf("\n the exit enthalpy is %f kj/kg or %f btu/lb",hout_s,hout_f);

               tdbo = hout_s/1.006;
               printf("\n tdbo is %f",tdbo);
           //    tdbo_f = hout_f/0.240;//deg f
               twbo = tdbo;


            svp_out = sat_vape(twbo);               //sat vap pressure calculated
            printf("\n sat vap p exit is %f Pa   svpn    %f",svp_out,svpn);



            sphu_out = sp_hum(svp_out,tdbo,twbo);

            printf("\n sp humidity exit is %f kg/kg of satureated air",sphu_out);


            printf("\n AT COOLING TOWER EXIT\n");


               spvol_out_s = sp_vol(tdbo,twbo,svp_out);//cum/kg
               spvol_out_f = spvol_out_s / 0.0624279606;//cuft/lb
               printf("\n the exit specific volume is %f cum/kg or %f cu ft/lb ",spvol_out_s,spvol_out_f);

               row_out_s = row(sphu_out,spvol_out_s);

               row_out_f = row_out_s/16.0184634;
               printf("\n the exit density is %f kg/cum or %f lb/cu ft",row_out_s,row_out_f);

               vel_in_s = (g_s * spvol_in_s)/aia ;//m/s
               vel_in_f = vel_in_s/0.00508;//fpm
               vel_out_s = (g_s*spvol_out_s)/cela ;
               vel_out_f =  vel_out_s/0.00508;

               row_fill_s = (row_in_s + row_out_s)/2 ;
               row_fill_f = (row_in_f + row_out_f)/2 ;
               spvol_fill_s = (spvol_in_s+spvol_out_s)/2;
               spvol_fill_f = (spvol_in_f+spvol_out_f)/2;   //

               vel_fill_s = (g_s*spvol_fill_s)/cela ;
               vel_fill_f = vel_fill_s/0.00508;

               printf("\n the specific volume at fill is %f cum/kg or %f cu ft/lb ",spvol_fill_s,spvol_fill_f);
               printf("\n the density at fill is %f kg/cum or %f lb/cu ft",row_fill_s,row_fill_f);
               printf("\n the velocity of air at inlet is %f m/s or %f fpm ",vel_in_s,vel_in_f);
               printf("\n the velocity of air at fill is %f m/s or %f fpm ",vel_fill_s,vel_fill_f);
               printf("\n the velocity of air at exit is %f m/s or %f fpm ",vel_out_s,vel_out_f);

               cfm_fan_s = g_s * spvol_out_s;//cum/s
               cfm_fan_f = g_f * spvol_out_f;
               printf("\n the volume flow rate of air at fan CFAN is %f cum/s or %f cu ft/min \n",cfm_fan_s,cfm_fan_s);



               pd_il_inh20 = (3*((2*pow(10,-7)*vel_in_f*vel_in_f)+(3*pow(10,-6)*vel_in_f))*(row_in_f/0.07));
        //       pd_il_pa = pd_il_inh20*249;//pascal
               printf("\n pd il %f in ",pd_il_inh20);

               pd_de_inh20 = (5* ((3*pow(10,-7)*vel_out_f*vel_out_f)-(9*pow(10,-6)*vel_out_f) + 0.003 )*(row_out_f/0.07));
        //       pd_de_pa = pd_de_inh20*249;//pascal
               printf("\n pd de %f in ",pd_de_inh20);



               pdf1 = (fh_f/fs_f);
               pdf2 =pow(pdf1,0.963077);
               pdf3 = ((7.44385*pow(10,-6))+((-4.078*pow(10,-6))*pow(fh_f,0.176686)));
               pdf4 = (4.897*(pow(10,-4))*(pow(vel_fill_f,2.3386)));
               pdf5 = (0.0428*wl_f);
               pdf6 = (pow(vel_fill_f,1.33558));

               pd_fill_inh20 = (pdf2*pdf3*(pdf4+(pdf5*pdf6)))*(row_fill_f/0.07);
               printf("\n pd fill %f in ",pd_fill_inh20);
         //      pd_fill_pa = pd_fill_inh20*249;

               tpd_s_inh20 = pd_il_inh20 + pd_de_inh20 + pd_fill_inh20;
               printf("\n hs %f in ",tpd_s_inh20);
               tpd_s_pa = tpd_s_inh20*249;
               pd_v_f = (row_out_f/1202312)*(cfm_fan_f/edisa_f)*(cfm_fan_f/edisa_f);
               printf("\n hv %f in",pd_v_f);
               pd_v_s = pd_v_f * 249 ;
               tpd_f = pd_v_f + tpd_s_inh20;
               tpd_s = tpd_s_pa + pd_v_s;
               printf("\n the overall pressure drop across the cooling tower is %f pa or %f inches of water col",tpd_s,tpd_f);

               bhp = ((cfm_fan_f*(tpd_s_inh20+pd_v_f))/(6356*fann));//hp
               bhp_w = bhp*0.745699872;//kw
               printf("\n the required brake horse power of the fan is %f Bhp or %f kw",bhp,bhp_w);

return (bhp_w);
}



float hooke(int max,float fh,float fs,float epsilon,float fhsl)
{

int iadj,iters,i,j,xs1,ys1,xs2,ys2,islfs = 1,fssl = 1;
float islfh,fhi,fhf,fsi,fsf,bhpi;

float a,b,m;
float bhp[250],fh1ta[250],fs1ta[250],bhp1ta[250],fh1pa,fs1pa,bhp1pa;
float fh1tn[250],fs1tn[250],bhp1tn[250],fh1pn,fs1pn,bhp1pn,fh1[250],fs1[250],bhp1[250],fh1_p,fs1_p,bhp1_p,fh1_m,fs1_m,bhp1_m;
float fh2ta[250],fs2ta[250],bhp2ta[250],fh2tn[250],fs2tn[250],bhp2tn[250];
float fh2pa,fs2pa,bhp2pa,fh2pn,fs2pn,bhp2pn,fh2_m,fs2_m,bhp2_m,fh2[250],fs2[250],bhp2[250];
float fhsp[250],fssp[250],lam1[250],lam2[250],lam[250],delta,fh3[250],fs3[250],bhp3[250],fh3_m,fs3_m,bhp3_m,bhp3_m0,bhp3_m1;


islfh = 0.001;
xs1 = 1;
ys1 = 0;
xs2 = 0;
ys2 = 1;
iadj = 0;
iters = 0;
fhi = 0.6096;
fhf = 1.8288;
fsi = 12.0f;
fsf = 19.0f;
a = (7.44385*pow(10,-6));
b = (-4.078*pow(10,-6));
m = 0.176686;


bhpi = bhpeqns(fhi , fsi);

fhi = fh;
fsi = fs;
bhpi = bhpeqns(fhi,fsi);
        i = 0;
        fh1ta[0] = fhi + (xs1 * islfh);
        fs1ta[0] = fsi + (ys1 * islfs);
        fh1pa = fh1ta[0] ;      fs1pa = fs1ta[0] ;
        bhp1pa = bhpeqns(fh1pa,fs1pa);
        bhp1ta[0] = bhp1pa ;
        fh1tn[0] = fhi - (xs1 * islfh);
        fs1tn[0] = fsi - (ys1 * islfs);
        fh1pn = fh1tn[0] ;      fs1pn = fs1tn[0] ;
        bhp1pn = bhpeqns(fh1pn,fs1pn);
        bhp1tn[0] = bhp1pn ;
        do
        {
            if(bhp1pa < bhpi)
            {
                j=0;
                while(bhp1pa < bhpi)
                {
                    fh1[0] = fhi + (xs1 * fhsl);
                    fs1[0] = fsi + ( ys1 * fssl);
                    fh1_p = fh1[0]  ;   fs1_p = fs1[0] ;
                    bhp1_p = bhpeqns(fh1_p,fs1_p);
                    bhp1[0] = bhp1_p;

                    if(bhp1_p < bhpi)
                    {
                       fh1[0] = fhi + (xs1 * fhsl);
                       fs1[0] = fsi + ( ys1 * fssl);
                       fh1_p = fh1[0] ;     fs1_p = fs1[0];
                     break;
                    }
                    else
                    {
                       fhsl = fhsl/2 ;
                       j++ ;
                    }
                }
                fh1[0] = fh1_p;
                fs1[0] = fs1_p;

            }
            else if(bhp1pn < bhpi)
            {
                j=0;
                do
                {
                    fh1[0] = fhi - (xs1 * fhsl);
                    fs1[0] = fsi - ( ys1 * fssl);
                    fh1_p = fh1[0]  ;   fs1_p = fs1[0] ;
                    bhp1_p = bhpeqns(fh1_p,fs1_p);
                    bhp1[0] = bhp1_p;

                    if(bhp1_p < bhpi)
                    {
                       fh1[0] = fhi - (xs1 * fhsl);
                       fs1[0] = fsi - ( ys1 * fssl);
                       fh1_p = fh1[0] ;     fs1_p = fs1[0];
                     break;
                    }
                    else
                    {
                       fhsl = fhsl/2 ;
                       j++ ;
                    }
                }while((bhp1pn < bhpi) && (fh1[0] >= fhi) && (fh1[0] <= fhf) && ( fs1[0] >= fsi) && (fs1[0] <= fsf));

                fh1[0] = fh1_p;
                fs1[0] = fs1_p;

            }
            else
            {
                fh1[0] = fhi ;
                fs1[0] = fsi ;
            }
        }while((fh1[0] >= fhi) && (fh1[0] <= fhf) && ( fs1[0] >= fsi) && (fs1[0] <= fsf));

        fh1_m = fh1[0] ;
        fs1_m = fs1[0] ;
        bhp1_m = bhpeqns(fh1_m,fs1_m);
        bhp1[0] = bhp1_m ;
        printf("\n  i = %d  fh1   %f   fs1   %f    bhp1    %f",i,fh1[0],fs1[0],bhp1[0]);

        fh2ta[0] = fh1_m + (xs2 * islfh);
        fs2ta[0] = fs1_m + (ys2 * islfs);
        fh2pa = fh2ta[0];   fs2pa = fs2ta[0];
        bhp2pa = bhpeqns(fh2pa , fs2pa);
        bhp2ta[0] = bhp2pa;

        fh2tn[0] = fh1_m - (xs2 * islfh);
        fs2tn[0] = fs1_m - (ys2 * islfs);
        fh2pn = fh2tn[0];   fs2pn = fs2tn[0];
        bhp2pn = bhpeqns(fh2pn , fs2pn);
        bhp2tn[0] = bhp2pn;

        if(bhp2pa < bhp1_m)
        {
            do
            {
                fh2[0] = fh1_m + (xs2 * fhsl);
                fs2[0] = fs1_m + (ys2 * fssl);
            }while((fh2[0] >= fhi) && (fh2[0] <= fhf) && ( fs2[0] >= fsi) && (fs2[0] <= fsf));
        }
        else if (bhp2pn < bhp1_m)
        {
            do
            {
                fh2[0] = fh1_m - (xs2 * fhsl);
                fs2[0] = fs1_m - (ys2 * fssl);
            }while((fh2[0] >= fhi) && (fh2[0] <= fhf) && ( fs2[0] >= fsi) && (fs2[0] <= fsf));
        }
        else
        {
            fh2[0] = fh1_m;
            fs2[0] = fs1_m;
        }

        fh2_m = fh2[0] ;
        fs2_m = fs2[0] ;
        bhp2_m = bhpeqns(fh2_m,fs2_m);
        bhp2[0] = bhp2_m ;
        printf("\n  i = %d  fh2   %f   fs2   %f    bhp2    %f",i,fh2[0],fs2[0],bhp2[0]);

        fhsp[0] = fh2_m - fh1_m;
        fssp[0] = fs2_m - fs1_m;
        lam1[0] = (-1*a)/(b*fh2_m*m);
        lam2[0] = 1/(m-1);
        lam[0] = fabs((pow(lam1[0],lam2[0]))-1) ;       //note absolute values only is added here
        delta = lam[0];
        printf("\n the value of direction search i  %d   delta   %f",i,delta);

        fh3_m = fh2_m + (lam[0] * fhsp[0]);
        fh3[0] = fh3_m;
        fs3_m = fs2_m + (1 * fssp[0]);
        fs3[0] = fs3_m;
        bhp3_m1 = bhpeqns(fh3_m,fs3_m);     ///*****
        bhp3_m0 = bhpeqns(fh3_m,fs2_m);
        if((bhp3_m1 < bhp3_m0) && (bhp3_m1 < bhp2_m))
        {
            bhp3_m = bhp3_m1;
            fh3[0] = fh3_m;
            fs3[0] = fs3_m;
        }
        else if((bhp3_m0 < bhp3_m1) && (bhp3_m0 < bhp2_m))
        {
            fh3[0] = fh3_m;
            fs3[0] = fs2_m;
            bhp3_m = bhp3_m0;
        }
        else
        {
            fh3[0] = fh2[0];
            fs3[0] = fs2[0];
        }
        bhp3[0] = bhp3_m;
        bhp[0] = bhp3[0];   fh = fh3[0];    fs = fs3[0];
        printf("\n  i = %d  fh3   %f   fs3   %f    bhp3    %f",i,fh3[i],fs3[i],bhp3[i]);


//for(i = 0 ; (i < max) && (bhp[i+1] <= bhp[i]) ; i++)
    {





i = 1;

do
{
    iadj++;
while(i<max)
{
iters++;
    fhi = fh3[i-1] ;
    fsi = fs3[i-1];


        fh1ta[i] = fhi + (xs1 * islfh);
        fs1ta[i] = fsi + (ys1 * islfs);
        fh1pa = fh1ta[i] ;      fs1pa = fs1ta[i] ;
        bhp1pa = bhpeqns(fh1pa,fs1pa);
        bhp1ta[i] = bhp1pa ;

        fh1tn[i] = fhi - (xs1 * islfh);
        fs1tn[i] = fsi - (ys1 * islfs);
        fh1pn = fh1tn[i] ;      fs1pn = fs1tn[i] ;
        bhp1pn = bhpeqns(fh1pn,fs1pn);
        bhp1tn[i] = bhp1pn ;

        do
        {
            if(bhp1pa < bhpi)
            {
                j=0;
                while(bhp1pa < bhpi)
                {
                    fh1[i] = fhi + (xs1 * fhsl);
                    fs1[i] = fsi + ( ys1 * fssl);
                    fh1_p = fh1[i]  ;   fs1_p = fs1[i] ;
                    bhp1_p = bhpeqns(fh1_p,fs1_p);
                    bhp1[i] = bhp1_p;

                    if(bhp1_p < bhpi)
                    {
                       fh1[i] = fhi + (xs1 * fhsl);
                       fs1[i] = fsi + ( ys1 * fssl);
                       fh1_p = fh1[i] ;     fs1_p = fs1[i];
                     break;
                    }
                    else
                    {
                       fhsl = fhsl/2 ;
                       j++ ;
                    }
                }
                fh1[i] = fh1_p;
                fs1[i] = fs1_p;

            }
            else if(bhp1pn < bhpi)
            {
                j=0;
                do
                {
                    fh1[i] = fhi - (xs1 * fhsl);
                    fs1[i] = fsi - ( ys1 * fssl);
                    fh1_p = fh1[i]  ;   fs1_p = fs1[i] ;
                    bhp1_p = bhpeqns(fh1_p,fs1_p);
                    bhp1[i] = bhp1_p;

                    if(bhp1_p < bhpi)
                    {
                       fh1[i] = fhi - (xs1 * fhsl);
                       fs1[i] = fsi - ( ys1 * fssl);
                       fh1_p = fh1[i] ;     fs1_p = fs1[i];
                     break;
                    }
                    else
                    {
                       fhsl = fhsl/2 ;
                       j++ ;
                    }
                }while((bhp1pn < bhpi) && (fh1[i] >= fhi) && (fh1[i] <= fhf) && ( fs1[i] >= fsi) && (fs1[i] <= fsf));

                fh1[i] = fh1_p;
                fs1[i] = fs1_p;

            }
            else
            {
                fh1[i] = fhi ;
                fs1[i] = fsi ;
            }
        }while((fh1[i] >= fhi) && (fh1[i] <= fhf) && ( fs1[i] >= fsi) && (fs1[i] <= fsf));

        fh1_m = fh1[i] ;
        fs1_m = fs1[i] ;
        fh1[i] = fhdet(fh1_m);
        fs1[i] = fsdet(fs1_m);
        fh1_m = fh1[i];
        fs1_m = fs1[i];

        bhp1_m = bhpeqns(fh1_m,fs1_m);
        bhp1[i] = bhp1_m ;

        printf("\n  i = %d  fh1   %f   fs1   %f    bhp1    %f",i,fh1[i],fs1[i],bhp1[i]);

        fh2ta[i] = fh1_m + (xs2 * islfh);
        fs2ta[i] = fs1_m + (ys2 * islfs);
        fh2pa = fh2ta[i];   fs2pa = fs2ta[i];
        bhp2pa = bhpeqns(fh2pa , fs2pa);
        bhp2ta[i] = bhp2pa;

        fh2tn[i] = fh1_m - (xs2 * islfh);
        fs2tn[i] = fs1_m - (ys2 * islfs);
        fh2pn = fh2tn[i];   fs2pn = fs2tn[i];
        bhp2pn = bhpeqns(fh2pn , fs2pn);
        bhp2tn[i] = bhp2pn;

        if(bhp2pa < bhp1_m)
        {
            do
            {
                fh2[i] = fh1[i] + (xs2 * fhsl);
                fs2[i] = fs1[i] + (ys2 * fssl);
            }while((fh2[i] >= fhi) && (fh2[i] <= fhf) && ( fs2[i] >= fsi) && (fs2[i] <= fsf));
        }
        else if (bhp2pn < bhp1_m)
        {
            do
            {
                fh2[i] = fh1[i] - (xs2 * fhsl);
                fs2[i] = fs1[i] - (ys2 * fssl);
            }while((fh2[i] >= fhi) && (fh2[i] <= fhf) && ( fs2[i] >= fsi) && (fs2[i] <= fsf));
        }
        else
        {
            fh2[i] = fh1_m;
            fs2[i] = fs1_m;
        }

        fh2_m = fh2[i] ;
        fs2_m = fs2[i] ;
        fh2[i] = fhdet(fh2_m);
        fs2[i] = fsdet(fs2_m);
        fh2_m = fh2[i];
        fs2_m = fs2[i];
        bhp2_m = bhpeqns(fh2_m,fs2_m);
        bhp2[i] = bhp2_m = bhp[i];

        printf("\n  i = %d  fh2   %f   fs2   %f    bhp2    %f",i,fh2[i],fs2[i],bhp2[i]);    

        fhsp[i] = fh2_m - fh1_m;
        fssp[i] = fs2_m - fs1_m;
        lam1[i] = (-1*a)/(b*fh2_m*m);
        lam2[i] = 1/(m-1);
        lam[i] = fabs((pow(lam1[i],lam2[i]))-1);       //note absolute values only is added here
        delta = lam[i];
        printf("\n the value of direction search i  %d   delta   %f",i,delta);

        fh3_m = fh2_m + (lam[i] * fhsp[i]);
        fh3[i] = fh3_m;
        fs3_m = fs2_m + (1 * fssp[i]);
        fs3[i] = fs3_m;
        bhp3_m1 = bhpeqns(fh3_m,fs3_m);     ///*****
        bhp3_m0 = bhpeqns(fh3_m,fs2_m);
        if((bhp3_m1 < bhp3_m0) && (bhp3_m1 < bhp2_m))
        {
            bhp3_m = bhp3_m1;
            fh3[i] = fh3_m;
            fs3[i] = fs3_m;
        }
        else if((bhp3_m0 < bhp3_m1) && (bhp3_m0 < bhp2_m))
        {
            fh3[i] = fh3_m;
            fs3[i] = fs2_m;
            bhp3_m = bhp3_m0;
        }
        else
        {
            fh3[i] = fh2[i];
            fs3[i] = fs2[i];
        }

        fh3_m = fh3[i];
        fs3_m = fs3[i];
        fh3[i] = fhdet(fh3_m);
        fs3[i] = fsdet(fs3_m);
        fh3_m = fh3[i];
        fs3_m = fs3[i];
        bhp3_m = bhpeqns(fh3_m,fs3_m);
        bhp3[i] = bhp3_m;
        bhp[i] = bhp3[i];
        fh = fh3[i];
        fs = fs3[i];

        printf("\n  i = %d  fh3   %f   fs3   %f    bhp3    %f",i,fh3[i],fs3[i],bhp3[i]);

        if (bhp[i] < bhp[i-1])
        {
            i++;
        }
        else
        {
            printf("\n the optimum value is reached at fh  %f    fs   %f    bhp    %f",fh,fs,bhp[i]);
        }
}
}while( (fh >= fhi) && (fh<= fhf) && (fs >= fsi) && (fs <= fsf) && (fhsl >= epsilon));
printf("\n the number of times the while loop is run is %d",iadj);

}
    for(i = 0 ; (i < imax) && (bhp[i+1] <= bhp[i]) ; i++)
    {

        printf("\n  i = %d  fh1   %f   fs1   %f    bhp1    %f",i,fh1[i],fs1[i],bhp1[i]);

        printf("\n  i = %d  fh2   %f   fs2   %f    bhp2    %f",i,fh2[i],fs2[i],bhp2[i]);
        printf("\n  i = %d  fh3   %f   fs3   %f    bhp3    %f",i,fh3[i],fs3[i],bhp3[i]);
    }
    printf("\n the optimized value is bhp   %f   fh %f    fs    %f ",bhp3_m,fh3[i],fs3[i]);
getch();
return (iters);
}

float fsdet(float fs)
{
    float fsi = 12.0f;
    float fsf = 19.0f; float fse;
    if(fs > fsf)
        {   fse = fsf;  }
        else
        {   fse = fs;   }
        if(fs < fsi)
        {   fse = fsi;  }
        else
        {   fse = fs;   }
return (fse);
}

float fhdet(float fh)
{
    float fhi = 0.6096;
    float fhf = 1.8288; float fhe;
    if(fh > fhf)
        {   fhe = fhf;  }
        else
        {   fhe = fh;   }
        if(fh < fsi)
        {   fhe = fhi;  }
        else
        {   fhe = fh;   }
return (fhe);
}

float sat_vap(float twbi)
{
    float too,sat[7],fsvp;

    too = twbi + to;

    sat[0] = svpp[0]/too;
    sat[1] = svpp[1] * log(too);
    sat[2] = svpp[3] * (twbi-0.01);
    sat[3] = svpp[2] * (pow(10,sat[2]));//
    sat[4] = svpp[5] * (1-(t9/too));
    sat[5] = svpp[4] * (pow(10,sat[4]));//
    sat[6] = sat[0] + sat[1] +sat[3] + sat[5] + svpp[6] ;
    sat[7] = pow(10,sat[6]);                //sat vap pressure calculated
    fsvp = sat[7];

return (fsvp);
}

float sp_hum(float svp,float tdbi,float twbi)
{

float hum[5],fsphum;

 hum[0] = sphup[0]/sphup[1];//
 hum[1] = sphup[2] * sphup[3] * (tdbi-twbi);
 hum[2] = svp - hum[1];
 hum[3] = (hum[0] * hum[2])/(sphup[2] - hum[2]); // sp humidity is calculated
 fsphum = hum[3];
 fsphum = fabs(fsphum);

return (fsphum);
}
float sat_vape(float twbo)
{
    float too,sat[7],fsvpe;

    too = twbo + to;

    sat[0] = svpp[0]/too;
    sat[1] = svpp[1] * log(too);
    sat[2] = svpp[3] * (twbi-0.01);
    sat[3] = svpp[2] * (pow(10,sat[2]));//
    sat[4] = svpp[5] * (1-(t9/too));
    sat[5] = svpp[4] * (pow(10,sat[4]));//
    sat[6] = sat[0] + sat[1] +sat[3] + sat[5] + svpp[6] + svpp[7] ;
    sat[7] = pow(10,sat[6]);                //sat vap pressure calculated
    fsvpe = sat[7];

return (fsvpe);
}

float sp_enthalpy(float tdbi,float sphum)
{
    float spe[5],fspe;
    spe[0] = spep[0] * tdbi ;
    spe[1] = (spep[1] * tdbi) + spep[2];
    spe[3] = sphum * spe[1] ;
    spe[4] = spe[0] + spe[3];
    fspe = spe[4];
    fspe = fabs(fspe);
return(fspe);
}
float sp_vol(float tdbi,float twbi,float svp)
{
    float vol[5],fspvol;
    vol[0] = spvolro * (tdbi + to);
    vol[1] = sphup[2] * sphup[3] * (tdbi - twbi);
    vol[2] = sphup[2] - svp + vol[1];
    vol[3] = sphup[1] * vol[2] ;
    vol[4] = (vol[0] * 1000)/vol[3] ;
    fspvol = vol[4];
    fspvol = fabs(fspvol);
return(fspvol);
}

float row(float sphum,float spvol)
{
    float row;
    row = (1 + sphum)/spvol;
    row = fabs(row);
return(row);
}

kindly help me with the code

sai.aleenus 0 Newbie Poster · Answer 18 · 2013-04-12T18:30:35+00:00

the initial thw and other such parameters are constant for a program and the varying values are celn,fh,fs

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 19 · 2013-04-12T21:00:45+00:00

kindly help me with the code

Difficult, given that the code as posted doesn't compile (at least not under GCC). I'll do what I can, but you can expect that a lot of changes will be needed.

Are you certain that the bhpeqns() function is correct? I'd hate to be fixing the wrong thing...

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 20 · 2013-04-12T22:13:14+00:00

I've checked it out, and the problem (or at least one problem) is indeed in bhpeqns(). On the second call to the function, it goes into an infinite loop in the first do/while() loop. Apparently, the input values for bhpeqns() on the second pass are out of bounds.

BTW, you are a) still using void main(), despite the fact that I said it was incorrect, and b) still using short variable names that don't help describe their purpose. This isn't Fortran 77, you aren't limited to 6 character variable names! Even Turbo C++ let's you use 32 character variables, so there is no good reason to limit yourself to these pitiful, cryptic abbreviations.

Oh, and there's still no reason to have the vapor pressure, humidity and enthalpy calculations inside of the first do/while() loop, as there is nothing in the loop that affects them. You can and should extract that whole section of code and move it to before the beginning of that loop. Even if you insist on using separate variables rather than arrays for them, you can at least do that much.

sai.aleenus 0 Newbie Poster · Answer 21 · 2013-04-13T06:42:48+00:00

ok...the bhpeqns() works fine when run individually but might be posing problem when called in hooke function

Schol-R-LEA 1,446 Commie Mutant Traitor Featured Poster · Answer 22 · 2013-04-13T13:21:09+00:00

I think what is happening is this: in the hooke() function as written, the first call to bhpeqns() uses the default parameters of (0.6096, 12.0). The second call then uses the parameters input by the user, which in this case are (1.0, 12.0). This second set of values seems to cause bhpeqns() to blow up. Since 1.0 is more or less in the middle of the range for the fh values, this should not be the case.

I'll keep investigating this as best as I can, but I'll be out of town most of today and tomorrow; I probably won't be able to get back to you until the day after tomorrow. Sorry.

sai.aleenus 0 Newbie Poster · Answer 23 · 2013-04-13T16:12:10+00:00

its ok....i would be really greatful to you if you could help me with this problem...by the way i really thank you for your valuable suggestions till date

sai.aleenus 0 Newbie Poster · Answer 24 · 2013-04-15T12:27:03+00:00

thank you very much.....
i tried another method for optimisation --- kuhn tucker conditions
i wrote a program using these conditions as when you pointed out that during the second loop run the prog code is blowing up... i checked and came to know that it might be due to linear nature of constraints..
i was able to optimise it but i am facing a very minor problem in one of the loop which i am posting here...please help with this code

void kkt_optimize()
{
    float ofh,ofs,obhp,bhp_max,bhp_min,fill_h[8],fill_s[8],bh_p[8],pd_fill[8],u1[4],u2[4],u3[4],u4[4],u5[4],u6[4],u7[4],u8[4],mu[4];
    float fhi = 0.6096,fsi = 12.0f,fhf = 1.8288, fsf = 19.0f;
    float cela,gpmpc,wl,wl_f;
    float fh0,fh1,fh2,fh3,fh4,fh5,fh6,fh7,fs0,fs1,fs2,fs3,fs4,fs5,fs6,fs7,vel0,vel1,vel2,vel3;
    float xd,xdd,yd,ydd,yddd,yde,bhp_0,bhp_1,bhp_2,bhp_3,bhp_4,bhp_5,bhp_6,bhp_7,bh_pp[10];
    int x0[10],xe[10],location = 0;
    float y0[10],ye[10];
    int gd = DETECT, gm;
     initgraph(&gd,&gm,"c:\\turbo\\tc\\bgi"); 

    cela = cell*celw ;
    gpmpc = gpm/celn ;
//  wl = gpmpc/cela;
    wl_f = (gpmpc*4.402867)/(cela*3.2808*3.2808);

    printf("\n condition 0: fh >= fhi");
    fill_h[0] = fhi;
    fh0 = fill_h[0];
    xd = (1 - pd_fill_n);
    xdd = 1/xd;
    fill_s[0] = pow(10,xdd);
    fs0 = fill_s[0];
    if((fill_s[0] >= fsi) && (fill_s[0] <= fsf))
    {
        fs0 = fill_s[0];
        printf("\n might be a possible soln ---- check for mu");
        //mu condition verfication
        vel0 = fill_velocity(fh0,fs0);
        u1[0] = (pd_fill_c * (pow(vel0,pd_fill_o)));
        u2[0] = (pd_fill_d * pow(wl_f,pd_fill_q));
        u3[0] = (pow(vel0,pd_fill_p));
        u4[0] = u1[0]*u2[0]*u3[0];// const
        u5[0] = pow(fs0,(-1*pd_fill_n));//
        u6[0] = (pd_fill_n * pd_fill_a * pow(fh0,(pd_fill_n - 1))); //+
        u7[0] = (pd_fill_n + pd_fill_m);
        u8[0] = (u7[0] * pd_fill_b * (pow(fh0, (u7[0] - 1))));//
        mu[0] = (-1)* (u4[0] * u5[0] * (u6[0] + u8[0]));

            if(mu[0] < 0 )
            {
                printf(" \n this condition solution is invalid as a constraint is violated");
                pd_fill[0] = 25;
                bh_p[0] = 50;
            }
            else
            {
                printf("\n this is a possible condition");
                pd_fill[0] = total_pd_eqns(fh0,fs0);
                bh_p[0] = bhpeqns(fh0,fs0);
            }
            bhp_0 = bh_p[0];
    }
    else
    {
        printf("constraint violated");
        fill_s[0] = 20;
        pd_fill[0] = 25;
        bh_p[0] = 50;
    }
    bhp_0 = bh_p[0];
}

which i am posting here...please help with this code wr which i am posting here...please help with this code which i am posting here...please help with this code whe which i am posting here...please help with this code which i am posting here...please help with this code wr which i am posting here...please help with this codee which i am posting here...please help with this code which i am posting here...please help with this code

the problem i am facing is with the output value of fill_s[0].....

sai.aleenus 0 Newbie Poster · Answer 25 · 2013-04-15T14:07:07+00:00

note:

define pd_fill_n 0.9673

define pd_fill_a -0.12932

define pd_fill_b 0.08896

define pd_fill_m 0.1055

define pd_fill_c -0.005738

define pd_fill_o 1.022

define pd_fill_d 0.0001775

define pd_fill_p 0.7706

define pd_fill_q 0.508

sai.aleenus 0 Newbie Poster · Answer 26 · 2013-04-15T16:03:43+00:00

as given in my previous post pd_fill_n is const...i forgot to add that part of the code previously so posted that after the 1st post...
pd_fill_n has been initialised properly and is working for the rest of the calc but encountering an error for fill_s[0]

hookes pattern search optimisaton

Recommended Answers Collapse Answers

All 32 Replies

Recommended Answers