topic Re: LP in SAS/IML Software and Matrix Computations

LP

ayaafarag — Wed, 09 Sep 2015 17:02:02 GMT

I have two questions concerning this code:

1- how to solve this problem using linear programming?

2- how to make these constraints equality constraints and how to add another 2 inequality constraints?

start fun(DEV) global (n,p);

sumf=0;

do i = 1 to 6;

sumf = sumf + (DEV[I]);

end;

return (sumf);

finish fun;

start con(DEV) global (lambda,lenght,thrd,estvar1,g,mseco,mspeco,i2,R,nx,n,ny,rxy,p);

c=j(3,1,0);

sumc1= DEV[1]-DEV[4]+(DEV[7]*NX[1,1])+(DEV[8]*NX[1,2]);

sumc2=DEV[2]-DEV[5]+(DEV[7]*NX[2,1])+(DEV[8]*NX[2,2]);

sumc3=DEV[3]-DEV[6]+(DEV[7]*NX[3,1])+(DEV[8]*NX[3,2]);

c[1]= NY[1]-sumc1;

c[2]=NY[2]-SUMC2;

c[3]=NY[3]-SUMC3;

return (c);

finish con;

call LP (rc, kres, "fun");

Re: LP

Rick_SAS — Wed, 09 Sep 2015 18:13:15 GMT

Firstly, do you have SAS/IML 13.1 or later? If so, you should consider using the LPSOLVE function, which is simpler to use.

From the documentation for the LP subroutine it appears you need to represent the problem a little differently than you have.

The objective function (in your case, the sum of the elements) is represented by a vector c and the objective function is c`*x.

The constraints are represented by a matrix A and a vector b, with the condition being A*x=b. (Although in the calling syntax, c is appended as a column to the A matrix.)

I'm not sure what you are trying to do, but here is some pseudo code that might help us communicate better.. I am guessing because you never defined nx or ny in your program.

proc iml;
nx = {1 2, 3 2, 4 0};   /* is this part of A? */
ny = {-1 0 1};          /* is this supposed to be b? */
c = j(8, 1, 1);          /* optimize c*x */
A = I(3) || -I(3) || nx; /* constraint */
b = ny;

Do these equations represent the problem you are trying to solve? Are there any inequalities in the constraints?

As I said, if you can run the LPSOLVE subroutine, you'll get to an answer much faster.

Re: LP

ayaafarag — Thu, 10 Sep 2015 09:25:10 GMT

I have defined nx and ny before writing the LP problem. I want to write this code inside proc iml command.

So what is the write command to run " LP" model in this case ?

And I need to put equality constraints but I do not know how to do this here?

here is my detailed code:

Proc iml;

option nonotes;

P=2 ; n=3 ; B= {1, 1};

X1= { 0.019831833

0.150408859

0.782582506

};

x2={ 1.008931805

0.365347241

0.055856515

};

Xt1 = t(x1);

xt2=t(x2);

xt=xt1 || xt2;

E= 1#normal(repeat(-1,n));

Y=xt* B+E;

my = sum(y)/n;

m1 = sum(xt1)/n;

m2 = sum(xt2)/n;

ny =(y-my)/sqrt(t(y-my)*(y-my));

nX1 =(xt1-m1)/sqrt (t(xt1-m1)*(xt1-m1));

nX2 =(xt2-m2)/sqrt(t(xt2-m2)*(xt2-m2));

nx= nx1||nx2;

start fun(DEV) global (n,p);

sumf=0;

do i = 1 to 6;

sumf = sumf + (DEV[I]);

end;

return (sumf);

finish fun;

start con(DEV) global (lambda,lenght,thrd,estvar1,g,mseco,mspeco,i2,R,nx,n,ny,rxy,p);

c=j(3,1,0);

sumc1= DEV[1]-DEV[4]+(DEV[7]*NX[1,1])+(DEV[8]*NX[1,2]);

sumc2=DEV[2]-DEV[5]+(DEV[7]*NX[2,1])+(DEV[8]*NX[2,2]);

sumc3=DEV[3]-DEV[6]+(DEV[7]*NX[3,1])+(DEV[8]*NX[3,2]);

c[1]= NY[1]-sumc1;

c[2]=NY[2]-SUMC2;

c[3]=NY[3]-SUMC3;

return (c);

finish con;

call LP (rc, kres, "fun");

Re: LP

ayaafarag — Thu, 10 Sep 2015 10:23:49 GMT

Am using SAS 9.2

Re: LP

ayaafarag — Thu, 10 Sep 2015 11:38:11 GMT

Here is my problem

x is 3*2 matrix

y is 3*1 vector

min dev(1)+dev(2)+dev(3)+dev(4)+dev(5)+dev(6)

s.t.

dev(7)*x11+dev(8)*x12+dev(1)-dev(4)=y1

dev(7)*x21+dev(8)*x22+dev(2)-dev(5)=y2

dev(7)*x31+dev(8)*x32+dev(3)-dev(6)=y3

dev(1),dev(2),dev(3),dev(4),dev(5),dev(6)>= 0

dev(7),dev(8) are unrestricted in sign

dev(1)*dev(4)=0

dev(2)*dev(5)=0

dev(3)*dev(6)=0

Re: LP

Rick_SAS — Thu, 10 Sep 2015 13:17:04 GMT

OK. Two comments:
1) there have been 7 releases of SAS/IML since 9.2, so you are missing out on a lot of cool stuff.
2) Your last three constraints are not linear, so this is not an LP problem. If those constraints are important, you might be able to use a nonlinear constrained optimization solver. The NLPNMS (Nelder-Mead Simplex Method) and NLPQN (Quasi-Newton Method) routines support nonlinear constraints and work very well.

Re: LP

ayaafarag — Thu, 10 Sep 2015 23:35:06 GMT

I tried the attached code, but I have two questions:

1-How to specify that the first three constraints and the last three are equality constraints and the rest are inequality constraints?

2-DO opt[10] and opt[11] are written correctly according to my problem?

Re: LP

ayaafarag — Sat, 12 Sep 2015 12:03:12 GMT

An extra question please the output gives me negative values for dev[1] and dev[5] however I specify dev[1]-dev[6] are nonnegative decision variables, so how to solve for this problem?

Re: LP

Rick_SAS — Sun, 13 Sep 2015 20:10:40 GMT

Use the BLC= option to specify a matrix that encodes the linear constraints. The first row specifies the lower bounds. The second row specifies the upper bounds. Subsequent rows specify linear constraints with equalities of inequalities. See the documentation page "Parameter Constraints" for the details about how to specify the constrant matrix. If I understand your problem, you will have 5 rows. The third through fifth rows will encode your linear constraints that are >= 0.

Use the NLC= option to specify the nonlinear constraints. I think your nonlinear constraints are

start nlcon(DEV);
   nlc=j(3,1,0);
   nlc[1]=dev[1]*dev[4];
   nlc[2]=dev[2]*dev[5];
   nlc[3]=dev[3]*dev[6];
   return (c);
finish nlcon;

These are all equality constraints, so you will use opt[10]=3 and opt[11]=3.

Re: LP

RobPratt — Wed, 23 Sep 2015 14:41:30 GMT

The nonlinear constraints are not needed here. The objective function will naturally drive them to be satisfied. In particular, any solution that has dev[1] > 0 and dev[4] > 0 can be strictly improved by reducing both dev[1] and dev[4] by the same amount (namely, min(dev[1],dev[4])) until one of them becomes 0. So you can stick with the linear formulation. Here's how you can do it with PROC OPTMODEL in SAS/OR:

   var dev {1..8};
   for {i in 1..6} dev[i].lb = 0;
   min z = sum {i in 1..6} dev[i];
   con lincon {i in 1..3}:
      dev[7] * x[i,1] + dev[8] * x[i,2] + dev[i] - dev[i+3] = y[i];

   solve;
   print dev;

Re: LP

ayaafarag — Fri, 25 Sep 2015 19:56:04 GMT

How define x matrix and y vector inside proc optmodel?

Re: LP

RobPratt — Fri, 25 Sep 2015 20:28:04 GMT

You can use the NUM statement in PROC OPTMODEL to define x and y. For example, here is x:

num x {1..3, 1..2} = [
0.019831833, 1.008931805,
0.150408859, 0.365347241,
0.782582506, 0.055856515
];

You could instead read the data from a SAS data set or use a formula in PROC OPTMODEL.

Re: LP

ayaafarag — Fri, 25 Sep 2015 21:26:47 GMT

if my problem is:

x is 3*2 matrix

y is 3*1 vector

min dev(1)+dev(2)+dev(3)+dev(4)+dev(5)+dev(6)

s.t.

dev(7)*x11+dev(8)*x12+dev(1)-dev(4)=y1

dev(7)*x21+dev(8)*x22+dev(2)-dev(5)=y2

dev(7)*x31+dev(8)*x32+dev(3)-dev(6)=y3

dev(1),dev(2),dev(3),dev(4),dev(5),dev(6)>= 0

dev(7),dev(8) are unrestricted in sign

so what is missing in the following code?

Proc iml;

option nonotes;

P=2 ; n=3 ; nruns=50; B= {1, 1};

X1= { 2.148910293

-1.756195507

-0.562194351

};

x2={ 0.558425545

-0.0446071

-3.180084241

};

nx=x1 || x2;

Do j=1 to nruns;

E1= 1#normal(repeat(-1,n-1));

E2=3#NORMAL(REPEAT(-1,1));

E=E1//E2;

y=nx* B+E;

obj={1 1 1 1 1 1 0 0};

coef = I(3) || -I(3) || nx;

rhs=y;

call lp(rc,dev,dual,coef,rhs);

　end;

Re: LP

ayaafarag — Wed, 30 Sep 2015 15:03:32 GMT

Below is my code, I have 3 questions:

1- how to formulate it as a minimum problem?

2- how to formulate the constraints as equality constraints?

3- I want the first 6 decision variables to be greater than or equal zero and the rest 2 decision variables to be unrestricted in sign, so how to do that here?

Proc iml;

option nonotes;

P=2 ; n=3 ; nruns=50; B= {1, 1};

X1= { 2.148910293

-1.756195507

-0.562194351

};

x2={ 0.558425545

-0.0446071

-3.180084241

};

Xt1 = t(x1);

xt2=t(x2);

nx=xt1 || xt2;

Do j=1 to nruns;

E1= 1#normal(repeat(-1,n-1));

E2=3#NORMAL(REPEAT(-1,1));

E=E1//E2;

nY=nx* B+E;

　/* the problem data */

obj = {1 1 1 1 1 1 0 0};

coef=I(3) || -I(3) || nx;

bb = {0}//ny;

/* embed the objective function in the coefficient matrix */

a = obj // coef;

/* solve the problem */

call lp(rc, x, dual, a, bb);

end;

Print devres ;

Re: LP

ayaafarag — Wed, 30 Sep 2015 17:02:13 GMT

Please, how to make the first two decision variables unrestricted in sign variables?