Statistical programming, matrix languages, and more

LP

Reply
Occasional Contributor
Posts: 17

LP

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");

SAS Super FREQ
Posts: 3,415

Re: LP

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.

Occasional Contributor
Posts: 17

Re: LP

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");

Occasional Contributor
Posts: 17

Re: LP

Am using SAS 9.2

Occasional Contributor
Posts: 17

Re: LP

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

SAS Super FREQ
Posts: 3,415

Re: LP

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.
Occasional Contributor
Posts: 17

Re: LP

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?

  

Occasional Contributor
Posts: 17

Re: LP

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?

SAS Super FREQ
Posts: 3,415

Re: LP

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.

 

SAS Employee
Posts: 424

Re: LP

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;

 

 

 

Occasional Contributor
Posts: 17

Re: LP

How define x matrix and y vector inside proc optmodel?

SAS Employee
Posts: 424

Re: LP

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.

Occasional Contributor
Posts: 17

Re: LP

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

 

 

Occasional Contributor
Posts: 17

Re: LP

  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;

 

 

 

 

 

 

 

Occasional Contributor
Posts: 17

Re: LP

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  ;

 

Ask a Question
Discussion stats
  • 14 replies
  • 1267 views
  • 0 likes
  • 3 in conversation