Minimize a function f(x) over a given domain - CLP optmodel

[harmonic]

Edited 2:
Yes sorry I modified the answer to add some informations.

1. I understood why is ifeasible, this is because there is no p value associated to the v value if I put the constraints. Anyway, why in the Problem Summary Section there is "Non linear" and then the solver MILP is linear?

3. I used solve linearize, it sometimes doesn't work if I change data.

5. "My goal is to use ......" did you mean the new question 1?

Anyway is that possible that with a restrict domain of p values for each r and a value it will find a function next to zero? It seems too optimistic.

[RobPratt]

1. That behavior (nonlinear objective in Problem Summary and MILP solver in Solution Summary) is expected when you use SOLVE LINEARIZE.  The original objective is nonlinear because of ABS, but the automated linearization converts to a linear problem and then calls the MILP solver.

3. Regarding "sometimes doesn't work":

Do you get an error message?  Do you get wrong results?  What SAS version are you using?

5. Whether the expectation of getting an optimal objective value close to zero is too optimistic depends on the data.  The solver will find such a solution if one exists.  A good exercise would be to independently verify the feasibility and objective value of the resulting solution returned by the solver.

[harmonic]

Thank you,

1. Is there a way to put the contraints of the values v and p but keeping just the couple of values that matches with r and a? because if I put just the con for p values or v values then it doesn't find the match in the table vdata or pdata.

2. if you launch the code with the csv that I passed you, does it work?

[RobPratt]

I ran with your csv files and got the following log messages:

ERROR: The expression for f was not linearized.
ERROR: Unable to linearize. The black-box solver might be suitable for this problem.
NOTE: Unable to create problem instance due to previous errors.

This error arises from the fact that some p values are missing in PDATA.csv.  For example:

35691501	PF		6228
35691502	PD		6227
35691503	PI		6229

You can avoid the error by modifying PDATA.csv or the PData data set or by using a WHERE= clause in the READ DATA statement:

   read data PData(where=(p ne .)) into REMI_CLASS_AOP=[remi class aop] p;

But then the resulting problem is infeasible because some <r,c> pairs have no corresponding <r,c,a> triple and the OneValue[r,c] constraint becomes 0 = 1.  One way to avoid this infeasibility is to update the REMI_CLASS set after reading PData:

   REMI_CLASS = REMI_CLASS inter (setof {<r,c,a> in REMI_CLASS_AOP} <r,c>);

After these changes, the SOLVE LINEARIZE statement runs successfully.

[harmonic]

I got the same issue Infeasible when I use this con because there is no corresponding value in the pData:

con Pcon {<r,c,a> in REMI_CLASS_AOP: not(9 <= p[r,c,a] < 11.5)}: X[r,c,a] = 0;

Is there a way to put constraints in v and p values and complete the algorithm successful?

[RobPratt]

If you don't want to modify the input data before calling PROC OPTMODEL, the next best approach is get rid of the disallowed <r,c,a> tuples as soon as possible by not even reading them:

   read data PData(where=(9 <= p < 11.5)) into REMI_CLASS_AOP=[remi class aop] p;

Now you don't need an explicit constraint to force the disallowed variables to 0 because those <r,c,a> tuples no longer even appear in the index set.

To make the resulting problem feasible, you can modify the REMI_CLASS set as before:

   REMI_CLASS = REMI_CLASS inter (setof {<r,c,a> in REMI_CLASS_AOP} <r,c>);

[harmonic]

My question is if I can use the "con" statement after reading the dataset. 

Because for example I want to add other variable to the vector REMI_CLASS_AOP, for example :

	set <str,str,str,num,num> REMI_CLASS_AOP;
	num p {REMI_CLASS_AOP};

    read data sas.PData into REMI_CLASS_AOP=[remi class aop q1 q2] p;

So I want to put a constraint like (q1-q2)/q2 < 0.05, but with the "con" statement, and not in the input dataset. I know that this doesn't make any sense but I would like to be sure that I can implement it. Perhaps the problem will be in the min function.

[RobPratt]

If q1 and q2 depend on <r,c,a>, I recommend the following:

	set <str,str,str> REMI_CLASS_AOP;
	num p {REMI_CLASS_AOP};
	num q1 {REMI_CLASS_AOP};
	num q2 {REMI_CLASS_AOP};
    read data sas.PData(where=((q1-q2)/q2 < 0.05)) into REMI_CLASS_AOP=[remi class aop] p q1 q2;

[harmonic]

Can I put it in a "constraint" statement?  like

con Pcon {<r,c,a> in REMI_CLASS_AOP: not(9 <= p[r,c,a] < 11.5)}: X[r,c,a] = 0;
con Qcon {<r,c,a> in REMI_CLASS_AOP: not((q1[r,c,a]-q2[r,c,a])/q2[r,c,a])<0.05}:
      X[r,c,a] = 0;

And don't have the infeasible problem? And keep just the couple of values that are in both v and p.

[RobPratt]

Yes, you can do that, but you need to move the final right parenthesis to correctly enforce the logical condition you want:

not((q1[r,c,a]-q2[r,c,a])/q2[r,c,a]<0.05)

To avoid infeasibility, you would still need to modify REMI_CLASS:

   REMI_CLASS = REMI_CLASS inter setof {<r,c,a> in REMI_CLASS_AOP: 9 <= p[r,c,a] < 11.5 and (q1[r,c,a]-q2[r,c,a])/q2[r,c,a]<0.05} <r,c>;

The alternative approach of omitting the Pcon and Qcon constraints and instead using the WHERE= option in the READ DATA statement would avoid the dual maintenance of checking these conditions in two different places.

[harmonic]

1. Where should I put this statement? 

 REMI_CLASS = REMI_CLASS inter setof {<r,c,a> in REMI_CLASS_AOP: 9 <= p[r,c,a] < 11.5 and (q1[r,c,a]-q2[r,c,a])/q2[r,c,a]<0.05} <r,c>;

Because it doen't make sense if the check would be done before of the "constraint" statement.

2. Can I put a constraint in the value of r as well? like r eq "30000301"?

thank you in advance

[RobPratt]

1. This statement should appear after the READ DATA PData and before the SOLVE statement.  The CON statements are just declarations and not programming statements, so it doesn't matter whether you update REMI_CLASS before or after the CON statements.

2. Yes, you can use the dummy index r as part of a logical condition in a constraint.

[harmonic]

1. What is wrong?

/* 117        solve linearize;
NOTE: Problem generation will use 4 threads.
ERROR: The array subscript 'v['34944801',SI]' is invalid at line 103 column 84.
NOTE: Unable to create problem instance due to previous errors.
118        
119        *print X;
120        *print {<r,c> in REMI_CLASS} (sum {<(r),(c),a> in REMI_CLASS_AOP} p[r,c,a]*X[r,c,a]);
121        
122        create data SolData1 from [r c a] v[r,c] p X;
ERROR: The array subscript 'v['34944801',SI]' is invalid at line 122 column 36. */


proc optmodel;
	set <str,str> REMI_CLASS; * dichiaro il tipo di variabile stringa;
	num v {REMI_CLASS}; * dichiaro la variabile numerica;

	read data sas.VData into REMI_CLASS=[remi class] v; * leggo il volume e lo inserisco in un vettore;

	set <str,str,str> REMI_CLASS_AOP;
	num p {REMI_CLASS_AOP};

    read data sas.PData into REMI_CLASS_AOP=[remi class aop] p;

	var X {REMI_CLASS_AOP} binary;
	
	REMI_CLASS = REMI_CLASS inter setof {<r,c,a> in REMI_CLASS_AOP: 9 <= p[r,c,a] < 11.5} <r,c>;
	
	/* funzione obiettivo */
	min f = abs(sum {<r,c,a> in REMI_CLASS_AOP} (if char(c,1) = 'M' then 1 else -1 ) *v[r,c]*p[r,c,a]*X[r,c,a]);
	
	/* constraints */
	con onevalue {<r,c> in REMI_CLASS}:
	   sum {<(r),(c),a> in REMI_CLASS_AOP} X[r,c,a] = 1; * Verrà posta a uno la variabile che minimizza;

	con Pcon {<r,c,a> in REMI_CLASS_AOP: 
		not(9 <= p[r,c,a] < 11.5)}: X[r,c,a] = 0;
	
	solve linearize;

	create data SolData1 from [r c a] v[r,c] p X;
	create data SolData2 from [r c a]={<r,c,a> in REMI_CLASS_AOP: X[r,c,a].sol > 0.5} v[r,c] p X;
	create data outdata_p from [remi class aop] p X ;
	create data outdata_v from [remi class] v ;
		
quit;

 
2. And if I would like to put constraint in v values like v < 100, how can I write the inter setof? Because it won't find the correspond in the pdata.

[RobPratt]

Both of those error messages mean that <'34944801',SI> is not in the index set for v but you are trying to access v['34944801',SI].  The first error message mentions line 103, which you did not show but should correspond to the MIN statement where v is used.  This error arises because the sum is over all <r,c,a> in REMI_CLASS_AOP, which you have not restricted to satisfy your condition on p.  To avoid this error, you can explicitly check membership in REMI_CLASS in the MIN statement:

   min f = abs(sum {<r,c,a> in REMI_CLASS_AOP: <r,c> in REMI_CLASS} (if char(c,1) = 'M' then 1 else -1)*v[r,c]*p[r,c,a]*X[r,c,a]);

To avoid the second error, you can make a similar change to the CREATE DATA statement:

   create data SolData1 from [r c a]={<r,c,a> in REMI_CLASS_AOP: <r,c> in REMI_CLASS} v[r,c] p X;

To restrict v < 100, you can modify the READ DATA statement:

   read data sas.VData(where=(v < 100)) into REMI_CLASS=[remi class] v;

The subsequent INTER SETOF does not then need modification.

[sbxkoenk]

Hello @harmonic ,

I de-duplicated your question by moving your topic thread in programming board to same topic thread in Operations Research board (this one!).

Instead of using SAS Optimization (definitely a good choice!), you can also do this via our Interactive Matrix Language (IML). From that IML module you also have access to several solvers.
However, I cannot do this short program quickly for you from the top of my head. Maybe @Rick_SAS can.

• Two tips for optimizing a function that has a restricted domain
  By Rick Wicklin on The DO Loop August 21, 2019
  https://blogs.sas.com/content/iml/2019/08/21/tips-optimize-restricted-domain.html
• Least-squares optimization and the Gauss-Newton method
  By Rick Wicklin on The DO Loop March 16, 2022
  https://blogs.sas.com/content/iml/2022/03/16/least-squares-gauss-newton.html

Koen