Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Options

- RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

☑ This topic is **solved**.
Need further help from the community? Please
sign in and ask a **new** question.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Posted 01-10-2023 03:43 PM
(461 views)

I am trying to write an optmodel algorithm using subproblems. I would like to solve for a parameter vector (lets call it V). For any V, I can estimate another vector W (using the nlp solver). W is a function of V -- W(V). I want the solution for V to minimize a function F=T-C*W(V). The problem I am having is that I cannot figure out how to make optmodel treat W(V) from the first problem as a function of V, and V itself does not explicitly appear in F=T-C*W(V). That is, SAS reads this as F=T-C*W with W just being the solution to the first problem not dependent on V.

1 ACCEPTED SOLUTION

Accepted Solutions

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Because you have the same variables and constraints at both levels, an alternative approach is to combine the two objectives and solve one QP. Explicitly, introduce a sufficiently large constant alpha, and minimize a weighted sum of your two objectives, where alpha is the weight on the first objective:

```
min Z = alpha * MSE + obj_W;
```

You would not need to use the PROBLEM, USE PROBLEM, or FIX statement. The idea is that the MSE part of the objective encourages W to take the "right" value with respect to V and obj_W is a secondary objective.

7 REPLIES 7

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Can you please show what you tried?

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Sure! Thank you for responding! I have pasted my optmodel code below with my comments.

proc optmodel Printlevel=0; /* Specifiy a variable with the number of control group rating areas */ number num_RA_control = %sysfunc(countw(&listRA_c.,%str( ))); number num_RA_all = num_RA_control+1; /* Sets to read in data */ set<num> NOBS ; /* number of observations (predictors) */ set<num> N_RAs = 1..num_RA_control ; /* Number of control rating areas */ set<num> N_RAs_tot = 1..num_RA_all ; /* Number of control rating areas + 1 (for the treated area) */ /* Row names */ string row_names {NOBS}; /* Standardized */ /* Control group matrix */ number control_matrix {NOBS,N_RAs}; read data std_matrix_notreat_B into NOBS=[_N_] row_names = _NAME_ %macro read_control; %do overRAs = 1 %to %sysfunc(countw(&listRA_c.,%str( ))); control_matrix[_N_,&overRAs.] = %scan(&listRA_c.,&overRAs.,%str( )) %end; ; %mend read_control; %read_control; /* print row_names control_matrix; */ /* Treated group vector */ number treated_vector {NOBS}; read data std_matrix_treat_B into NOBS=[_N_] treated_vector = /*%scan(&listRA_t.,&i.,%str( ))*/RA_FL01 ; /* print row_names treated_vector; */ /* Treated + Controls Matrix (for setting the initial guess of V) */ number matrix_all_bronze {NOBS,N_RAs_tot}; for {i in NOBS} matrix_all_bronze[i,1]=treated_vector[i]; for {i in NOBS, j in 2..num_RA_all} matrix_all_bronze[i,j]=control_matrix[i,j-1]; /* Create initial vectors for weights (W - m rating areas X 1) and the weighting matrix (V - n predictors x n predictors) */ /* initialize all weights to be equal and sum to 1 */ number init_vec_W {N_RAs} init 1/num_RA_control; /* std dev of the predictors */ number init_vec_V {NOBS}; for {i in NOBS} init_vec_V[i] = ((sum {j in N_RAs_tot} (matrix_all_bronze[i,j]-((sum {k in N_RAs_tot} matrix_all_bronze[i,k])/num_RA_all))**2)/(num_RA_all-1))**0.5; for {i in NOBS} init_vec_V[i] = init_vec_V[i]/init_vec_V[9]; /* For any given V we can calculate W - This is the easy part of the problem */ number Wsol {i in N_RAs,j in 1..3} init 0; /* keep track of different solutions for W */ /* two approaches */ /* 1. Use only the MSE (assumes V=1 for all v) */ var parm_W {i in N_RAs} init init_vec_W[i] >=0 <=1; con w_con1 : sum {i in N_RAs} parm_W[i] = 1; min MSE = sum {i in NOBS} ((treated_vector[i] - sum {j in N_RAs} control_matrix[i,j]*parm_W[j]))**2; problem v1 include parm_W MSE w_con1; use problem v1; solve with qp; for {i in N_RAs,j in 1..2} Wsol[i,1]=parm_W[i]; /* 2. Use the full objective function (assumes V=init_vec_V) */ number V {i in NOBS} = init_vec_V[i]; min obj_W = sum {i in NOBS} (treated_vector[i] - sum {j in N_RAs} control_matrix[i,j]*parm_W[j])*(V[i])*(treated_vector[i] - sum {k in N_RAs} control_matrix[i,k]*parm_W[k]); problem v2 include parm_W obj_W w_con1; solve with qp; for {i in N_RAs,j in 1..2} Wsol[i,2]=parm_W[i]; /* From the above problems you can see that W is a function of V as W(v1)~=W(v2) of the two approaches above v2 is better (obj_W<MSE) */ /* The above methods do not get me close to getting an estimate for V. To solve for V I need to minimize the the MSE objective function with respect to V. As an intermediate step, for any iteration of a V value we can calculate W (as done above) Then apply that W(V) to the MSE objective function. I want to find the V that minimizes the objective MSE where W(V) is calculated by minimizing obj_W. For example: In IML I can specify a function that I then call to minimize: start synth(v2) global(X1_SF,X0_SF,Y1_ESI_SF,Y0_ESI_SF,F1_ESI_SF,F0_ESI_SF,X_ns_SF,Y_ESI_ns_SF); D = 1//v2`; V = diag(D); ** Quadratic expression ** H = F0_ESI_SF`*V*F0_ESI_SF; ** Linear expression ** f = -F1_ESI_SF`*V*F0_ESI_SF; blc1 = j(1,ncol(F0_ESI_SF),0)||.||.; blc2 = j(1,ncol(F0_ESI_SF),1)||.||.; blc3 = j(1,ncol(F0_ESI_SF),1)||0||1; blc = blc1//blc2//blc3; init_0 = j(1,ncol(F0_ESI_SF),1/ncol(F0_ESI_SF)); optn1 = {0,0}; call NLPQUA(rc,W,H,init_0,optn1,blc,,,,f); e = F1_ESI_SF- F0_ESI_SF*W`; error = sum(e##2); return(error); AD = X_ns_SF//Y_ESI_ns_SF; s = sqrt(var(AD`)); v_start = s/s[1]; v_start = v_start[2:ncol(v_start)]; optn = {0,0}; con1 = j(1,nrow(v_start),0); con2 = j(1,nrow(v_start),.); con = con1//con2; call nlpnra(rc,result_ESI_SF,"synth",v_start,optn,con); finish; I am looking to do the exact same thing in OR. */ /* I tried writing this as a subproblem but I could not figure out how to make W a function of V */ /* My failed attempt below (because I just cannot figure out how to make W a function of V) NOTE: Iterating between the objective function obj_W where we hold V fixed then hold W fixed is not a way to solve this problem as V will be pushed to 0 for all values. */ num epsilon = 1E-6; num iter init 0; num exit_crit; num constant_V {NOBS}; var parm_V {i in NOBS} init init_vec_V[i] >=0; min obj_Wx = sum {i in NOBS} (treated_vector[i] - sum {j in N_RAs} control_matrix[i,j]*parm_W[j])*(parm_V[i])*(treated_vector[i] - sum {k in N_RAs} control_matrix[i,k]*parm_W[k]); min MSEx = sum {i in NOBS} ((treated_vector[i] - sum {j in N_RAs} control_matrix[i,j]*parm_W[j]))**2; problem Wx include parm_W parm_V obj_Wx w_con1; problem Vx include parm_W parm_V MSEx; /* main loop */ do while (1); /* Iterations */ iter = iter+1; /* Keep track of V */ for {i in NOBS} constant_V[i] = parm_V[i]; /* Solve for RA weights */ use problem Wx; fix parm_V; solve with nlp;/* optmodel does not think this is quadratic with V fixed */ /* Solve for V (how do I make the W from problem W a function of V for this problem?) */ use problem Vx; unfix parm_V; fix parm_W; solve with nlp; unfix parm_W; /* I would want this to calculate V - not W, but I cannot figure that out. */ /* Lets assume I got it to calculate V, then I would do this: */ /* Set an exit criteria such that when V is not changing much we exit (and calculate the final W) */ exit_crit = sum {i in NOBS} (parm_V[i]-constant_V[i])**2; print "Iteration:" iter "Exit criteria value:" exit_crit "Objective value W:" obj_Wx "Objective value MSE:" MSEx; if (exit_crit <= epsilon) or (iter>=50) then do; print "Estimation complete at iteration " iter; use problem Wx; fix parm_V; solve with nlp; print "Objective value:" obj_Wx; leave; end; end; for {i in N_RAs,j in 1..3} Wsol[i,3]=parm_W[i]; num synthetic_control {NOBS}; num treat_v_synth_dif {NOBS}; for {i in NOBS} synthetic_control[i] = (sum {j in N_RAs} control_matrix[i,j]*parm_W[j]); for {i in NOBS} treat_v_synth_dif[i] = (treated_vector[i] - sum {j in N_RAs} control_matrix[i,j]*parm_W[j]); print row_names treated_vector synthetic_control treat_v_synth_dif; /* parm_V and parm_W are our estimates - we can now apply the parm_W to all of the data to generate and outcome */ set<num> NOBS_allyrs; /* number of observations on control data */ number prems_all_years_control {NOBS_allyrs,N_RAs}; string row_names_allyrs {NOBS_allyrs}; read data prem_forOR_allyears_control_B into NOBS_allyrs=[_N_] row_names_allyrs = _NAME_ %macro read_control; %do overRAs = 1 %to %sysfunc(countw(&listRA_c.,%str( ))); prems_all_years_control[_N_,&overRAs.] = %scan(&listRA_c.,&overRAs.,%str( )) %end; ; %mend read_control; %read_control; number prems_all_years_treat {NOBS_allyrs}; read data prem_forOR_allyears_treat_B into NOBS_allyrs=[_N_] prems_all_years_treat = RA_FL01 ; num synthetic_control_prems {NOBS_allyrs}; num treat_v_synth_dif_prems {NOBS_allyrs}; for {i in NOBS_allyrs} synthetic_control_prems[i] = (sum {j in N_RAs} prems_all_years_control[i,j]*parm_W[j]); for {i in NOBS_allyrs} treat_v_synth_dif_prems[i] = (prems_all_years_treat[i] - sum {j in N_RAs} prems_all_years_control[i,j]*parm_W[j]); print row_names_allyrs prems_all_years_treat synthetic_control_prems treat_v_synth_dif_prems; quit;

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Thank you!

You are correct that this is a bilevel optimization problem and it seems that optmodel would not allow for me to solve that type of problem.

It does look like your suggestion would work in general for my problem. Unfortunately the data is confidential and lives on an air-gapped network so I would not be able to use SAS Optimization and connect to the cloud. It seems like I would need to go back to IML for this problem. Outside of IML, is there anything in SAS that allows one to solve bilevel optimization problems?

Thank you again!!

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Because you have the same variables and constraints at both levels, an alternative approach is to combine the two objectives and solve one QP. Explicitly, introduce a sufficiently large constant alpha, and minimize a weighted sum of your two objectives, where alpha is the weight on the first objective:

```
min Z = alpha * MSE + obj_W;
```

You would not need to use the PROBLEM, USE PROBLEM, or FIX statement. The idea is that the MSE part of the objective encourages W to take the "right" value with respect to V and obj_W is a secondary objective.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Thank you for all of your very quick replies!!

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

I'm always glad to help. An alternative approach to handle two objectives when you have a priority order is to solve twice, first with the primary objective and then with an "objective cut" and the secondary objective:

```
solve obj MyPrimaryObjective;
num best;
best = MyPrimaryObjective.sol;
con ObjectiveCut:
MyPrimaryObjective <= best; /* change to >= if maximization */
solve obj MySecondaryObjective;
```

Registration is open! SAS is returning to Vegas for an AI and analytics experience like no other! Whether you're an executive, manager, end user or SAS partner, SAS Innovate is designed for everyone on your team. Register for just $495 by 12/31/2023.

**If you are interested in speaking, there is still time to submit a session idea. More details are posted on the website. **

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.