Here is the whole model:
proc fcmp outlib=work.funcs.f;
function f1(Q, H);
return (
(7.18895803108255)+(Q)*(0.90898166129053)+(H)*(-0.52518338771685)+(Q**2)*(-0.07243939757267)+(Q)*(H)*(-0.0101157174044)+(H**2)*(0.01140027697845)+(Q**3)*(0.00271247731414)+(Q**2)*(H)*(0.00078575662965)+(Q)*(H**2)*(-8.0331130545237E-7)+(H**3)*(-0.00010011945215)+(Q**4)*(-0.00005707623607)+(Q**3)*(H)*(-0.0000145464761)+(Q**2)*(H**2)*(-4.0740624504878E-6)+(Q)*(H**3)*(5.6715539778245E-7)+(H**4)*(3.0399728809992E-7)+(Q**5)*(6.2706019750555E-7)+(Q**4)*(H)*(-4.9275589970883E-8)+(Q**3)*(H**2)*(1.1402170133176E-7)+(Q**2)*(H**3)*(-1.1450221400824E-8)
);
endsub;
function f2(Q, H);
return (
(9.56532591440209)+(Q)*(1.56477959347719)+(H)*(-0.58834848839521)+(Q**2)*(-0.14923333670627)+(Q)*(H)*(-0.02042582402532)+(H**2)*(0.01042030607175)+(Q**3)*(0.00551139580932)+(Q**2)*(H)*(0.00200710134061)+(Q)*(H**2)*(0.00006229314138)+(H**3)*(-0.00006493899805)+(Q**4)*(-0.00009774971067)+(Q**3)*(H)*(-0.00004756024059)+(Q**2)*(H**2)*(-0.00001062838926)+(Q)*(H**3)*(3.7950591182214E-7)+(Q**5)*(8.5938390272548E-7)+(Q**4)*(H)*(1.7658335207072E-7)+(Q**3)*(H**2)*(2.2050828738805E-7)+(H**5)*(8.7276517928309E-10)
);
endsub;
function g1(Q1,Q2); return ( 0.003134 * ((Q1+Q2)**2) + 0.003089 * (Q1**2) ); endsub;
function g2(Q1,Q2); return ( 0.003134 * ((Q1+Q2)**2) + 0.003089 * (Q2**2) ); endsub;
run;
%let Q1max = 30.4;
%let Q2max = 29.5;
%let P1min = 2.5;
%let P2min = 3.3;
%let P1max = 21.7;
%let P2max = 21.6;
%let PSETL = 1;
%let PSETR = 43;
%let HSETL = 70;
%let HSETR = 94;
%let STEP = 0.5;
options cmplib=work.funcs;
proc optmodel nthreads=64;
number P;
number H;
var Q1 >= 0 <= &Q1max;
var Q2 >= 0 <= &Q2max;
impvar HL1 = H - g1(Q1,Q2);
impvar HL2 = H - g2(Q1,Q2);
impvar P1 = Q1 * HL1 * f1(Q1, HL1) * 9.81 / 1000;
impvar P2 = Q2 * HL2 * f2(Q2, HL2) * 9.81 / 1000;
minimize Q = Q1 + Q2;
constraint P = P1 + P2;
constraint 0 <= P1 <= &P1max; /* (*) */
constraint 0 <= P2 <= &P2max; /* (*) */
set PSET = &PSETL .. &PSETR by &STEP;
set HSET = &HSETL .. &HSETR by &STEP;
num Qopt {PSET, HSET};
num Q1opt {PSET, HSET};
num Q2opt {PSET, HSET};
num P1opt {PSET, HSET};
num P2opt {PSET, HSET};
str solstatus {PSET, HSET};
num opterror {PSET, HSET};
do P = PSET;
do H = HSET;
solve with nlp / maxiter=10000 multistart=(maxstarts=10) seed=54321;
put P= H= HL1= HL2= P1= P2= Q1= Q2= _solution_status_=;
Qopt[P,H] = Q.sol;
Q1opt[P,H] = Q1.sol;
Q2opt[P,H] = Q2.sol;
P1opt[P,H] = P1;
P2opt[P,H] = P2;
solstatus[P,H] = _solution_status_;
opterror[P,H] = _OROPTMODEL_NUM_['OPTIMALITY_ERROR'];
end;
end;
create data Qdata from [P H] P1opt P2opt Q1opt Q2opt Qopt solstatus opterror;
quit;
Look at the rows marked with (*)
constraint 0 <= P1 <= &P1max; /* (*) */
constraint 0 <= P2 <= &P2max; /* (*) */
As you can see, P1 and P2 have 0 as lower bound, and it worked very well. Now I'm asked to restrict P1 and P2 interval from [0,PMAX] to [PMIN, PMAX], but zero must remain allowed, so the constraints should become
(P1 = 0) OR (&P1min <= P1 <= &P1max)
(P2 = 0) OR (&P2min <= P2 <= &P2max)
Hence this question. I hope I made myself clear now!
Thanks again,
Regards
... View more