$onText
Minimize the power loss during the operation of a hydrostatic thrust
bearing subject to a number of constraints.

References:
J.N. Siddall, Optimal Engineering Design. Marcel Dekker, New York, 1982.

K. Deb, M. Goyal, Optimizing engineering designs using a combined
genetic search. In: Back, T., (Ed.) Proceedings of the seventh International
Conference on Genetic Algorithms, 1997, East Lansing, MI, USA, July 19-23,
1997, pp.521-528.

Carlos A. Coello Coello, Treating constraints as objectives for single-
objective evolutionary optimization. Engineering Optimization, vol.32,
(2000), pp.275-308.
$offText

Scalars
  gamma            'weight density of oil (lb/in^3)'         /0.0307/
  C                'specific heat of oil (Btu/lb degr F)'    /0.5/
  N                'angular speed of shaft (rpm)'            /750/
  Ws               'thrust load (lb)'                        /101000/
  Pmax             'maximum inlet pressure (psi)'            /1000/
  delta_t_max      'maximum temperature change (degF)'       /50/
  hmin             'minimum film thickness (in)'             /0.001/
  g                'gravitational constant (in/s^2)'         /386.4/
  taf              'ambient temperature (degF)'              /100.0/
  P1               'outlet pressure, atmosferic (psi gauge)' /0/
  pump_efficiency  '70% efficiency of pump'                  /0.7/
  pi               'famous constant' ;

  pi = 2*arctan(inf);

Table oil_constants(*,*) oil grades

                    C1        n
    'SAE 5'       10.85     -3.91
    'SAE 10'      10.45     -3.72
    'SAE 20'      10.04     -3.55
    'SAE 30'       9.88     -3.48
    'SAE 40'       9.83     -3.46
    'SAE 50'       9.82     -3.44 ;

Parameters
  C1    'c factor for given oil grade'
  cn    'n factor for given oil grade'
  gr    'specific gravity of oil'
  tar   'ambient temperature'  ;

* Compute or extract parameters:
  C1 = oil_constants('SAE 20', 'C1');
  cn = oil_constants('SAE 20', 'n');
  gr = gamma/0.0361111;
  tar= 459.7 + taf;

Variables
  R         'bearing step radius (in)'
  R0        'recess radius (in)'
  mu        'oil viscosity in reynolds (Lb sec/in)'
  Q         'flow rate of oil (in^3/sec)'
  PL        'power loss (in lb/sec)'
  P0        'inlet pressure (psi gauge)'
  Ef        'power loss due to friction (in lb/sec)'
  Ep        'pumping energy (in lb/sec)'
  h         'film thickness (in)'
  delta_t   'temperature rise (deg F)'
  t         'temperature (deg F)'
  W         'weight (lb)' ;

Equations
  power_loss      'objective function'
  pumping_energy  'pump efficiency'
  friction        'friction loss'
  temp_rise       'temperature rise of the oil passing through the bearing'
  inlet_pressure  'relation between inlet pressuire and flow rate'
  load_capacity   'load-carrying capacity'
  oil_viscosity   'oil viscosity as a function of temperature'
  temperature     'T = ambient and oil outlet temperatures'
  radius          'inner radius should be the smallest'
  limit1          'limit exit loss to 0.1% of the pressure drop'
  limit2          'avoid surface damage in case of pressure loss' ;

* Objective function:
power_loss..  PL =e= Ep + Ef;

* Constraints:
pumping_energy.. Ep =e= Q*(P0-P1)/pump_efficiency;

friction..  Ef*h =e= sqr(2*pi*N/60)*[(2*pi*mu)]*(R**4-R0**4)/4;

temp_rise.. delta_t*(9336*Q*gamma*C) =e= Ef;

load_capacity.. W*(log(R)-log(R0)) =e= [(pi*P0)/2]*(sqr(R)-sqr(R0));

inlet_pressure.. P0*(pi*h**3) =e= (6*mu*Q)*(log(R)-log(R0));

oil_viscosity.. log10(8.112e6*mu+0.8) =e= (T**cn)*(10**C1);

temperature.. T =e= 560 + delta_t/2;

radius.. R =g= R0;

limit1.. gamma*sqr(Q) =l= 0.001*g*sqr(2*pi*R*h)*P0;

limit2.. W =l= 5000*pi*(sqr(R)-sqr(R0));

* Bounds:
P0.up = pmax;
h.lo = hmin;
delta_t.up = delta_t_max;
w.lo = Ws;

Ep.lo = 1;
Ef.lo = 1;
P0.lo = 1;
R.lo = 1;        R.up = 16;
R0.lo = 1;       R0.up = 16;
Q.lo = 1;        Q.up = 16;
mu.lo = 1.0e-6;  mu.up = 16.0e-6;
t.lo = 100;

* Variable scaling:
mu.scale = 1.0e-6;
h.scale = hmin;
W.scale = Ws;
PL.scale = 1.0e4;
Ep.scale = 1.0e4;
Ef.scale = 1.0e4;

* Initial point:
R.l = 6;
R0.l = 5;
mu.l = 0.000006;
Q.l = 3;
P0.l = 1000.0000;
Ef.l = 16000;
Ep.l = 3000;
h.l = 0.001;
delta_t.l = 50;
t.l = 600;
W.l = 101000;

Model htb /all/;
htb.scaleopt = 1;

Solve htb minimizing PL using nlp;

$ifThenI x%mode%==xbook
file out /htb.dat/
put out;
put R.l:15:7,   R0.l:15:7,  mu.l:15:7,  Q.l:15:7 /;
put P0.l:15:7,  Ef.l:15:7,  Ep.l:15:7,  h.l:15:7 /;
put delta_t.l:15:7,  t.l:15:7, W.l:15:7 /;
put /"Objective ="  PL.l:15:7 /;
$endIf
* End of htb