Optmodel algorithm

sml55 — Tue, 10 Jan 2023 20:43:37 GMT

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.

Re: Optmodel algorithm

RobPratt — Wed, 11 Jan 2023 03:05:27 GMT

Can you please show what you tried?

Re: Optmodel algorithm

sml55 — Wed, 11 Jan 2023 15:26:18 GMT

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;

Re: Optmodel algorithm

RobPratt — Wed, 11 Jan 2023 19:11:32 GMT

If you can write W as an explicit function of V, you can use an IMPVAR statement to declare W as an implicit variable. But your use case instead appears to be a bilevel optimization problem, which you can solve by using the CASEVALN function with the black-box optimization solver in the runOptmodel action in SAS Optimization. We have two documentation examples to illustrate this approach. The first one uses the glm action to solve the inner problem, and the second one uses the tsp action to solve the inner problem. In both cases, the SOLVE statement calls the black-box optimization solver from the runOptmodel action to solve the outer problem. For your use case, you can use either the solveQp action or the runOptmodel action to solve the inner problem.

Re: Optmodel algorithm

sml55 — Wed, 11 Jan 2023 20:22:05 GMT

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!!

Re: Optmodel algorithm

RobPratt — Wed, 11 Jan 2023 22:01:56 GMT

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.

Re: Optmodel algorithm

sml55 — Thu, 12 Jan 2023 13:55:10 GMT

Thank you for all of your very quick replies!!

Re: Optmodel algorithm

RobPratt — Thu, 12 Jan 2023 14:51:49 GMT

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;

topic Re: Optmodel algorithm in Mathematical Optimization, Discrete-Event Simulation, and OR

Optmodel algorithm

Re: Optmodel algorithm

Re: Optmodel algorithm

Re: Optmodel algorithm

Re: Optmodel algorithm

Re: Optmodel algorithm

Re: Optmodel algorithm

Re: Optmodel algorithm