If I understand what you are doing, you are using simulation from MVN to generate D instances of the random effects. Then you are using NLPNRR to fit a MLE model for each instance and are averaging the result.
That's not how mixed models are fit. As far as I know, your attempt to evaluate the "log-simulated likelihood" is not a valid estimation method. (And if you are re-generating the random effects between consecutive fits, this process probably won't converge.)
To answer your programming question, you can use the BLC option in the NLPNRR function to set the Boundary and Linear Constraints (BLC) for the parameters. In particular, you can specify that certain parameters must be positive. For an example of using BLC, see "Step 2: Set up constraints" in the article "Maximum likelihood estimation in SAS/IML."
For an example that compares the PROC IML approach to MLE with the PROC NLMIXED approach, see "Two ways to compute maximum likelihood estimates in SAS." I believe you can actually use PROC GLIMMIX for this problem, but I am not an expert on GLIMMIX and I am not 100% clear about what model you are trying to fit,
I suggest you read about "Types of logistic (or logit) models that can be fit using SAS" and focus on the models that include random effects. You might also want to browse the excellent papers by Robin High, in case they are relevant for you:
"Models for Ordinal Response Data"
"Fitting Statistical Models with PROCs NLMIXED and MCMC"
... View more