title: Optimization with Constraints
body: <p>I am working with the output from a model in which there are parameter estimates that may not follow a-priori expectations. I would like to write a function that forces these utility estimates back in line with those expectations. To do this, the function should minimize the sum of the squared deviance between the starting values and the new estimates. Since we have a-priori expections, the optimization should be subject to the following constraints:</p>
  
  <pre><code>B0 &lt; B1
  B1 &lt; B2
  ...
  Bj &lt; Bj+1 
  </code></pre>
  
  <p>For example, the raw parameter estimates below are flipflopped for B2 and B3. The columns <code>Delta</code> and <code>Delta^2</code> show the deviance between the original parameter estimate and the new coefficient. I am trying to minimize the column <code>Delta^2</code>. I've coded this up in Excel and shown how Excel's Solver would optimize this problem providing the set of constraints:</p>
  
  <pre><code>Beta    BetaRaw    Delta    Delta^2    BetaNew
  B0       1.2       0        0          1.2
  B1       1.3       0        0          1.3
  B2       1.6       -0.2     0.04       1.4
  B3       1.4       0        0          1.4
  B4       2.2       0        0          2.2
  </code></pre>
  
  <p>After reading through <code>?optim</code> and <code>?constrOptim</code>, I'm not able to grok how to set this up in R. I'm sure I'm just being a bit dense, but could use some pointers in the right direction!</p>
  
  <p>3/24/2012 - Added bounty since I'm not smart enough to translate the first answer.</p>
  
  <p>Here's some R code that should be on the right path. Assuming that the betas start with:</p>
  
  <pre><code>betas &lt;- c(1.2,1.3,1.6,1.4,2.2)
  </code></pre>
  
  <p>I want to minimize the following function such that <code>b0 &lt;= b1 &lt;= b2 &lt;= b3 &lt;= b4</code></p>
  
  <pre><code>f &lt;- function(x) {
    x1 &lt;- x[1]
    x2 &lt;- x[2]
    x3 &lt;- x[3]
    x4 &lt;- x[4]
    x5 &lt;- x[5]
  
   loss &lt;- (x1 - betas[1]) ^ 2 + 
           (x2 - betas[2]) ^ 2 + 
           (x3 - betas[3]) ^ 2 + 
           (x4 - betas[4]) ^ 2 +
           (x5 - betas[5]) ^ 2    
  
    return(loss)
  }
  </code></pre>
  
  <p>To show that the function works, the loss should be zero if we pass the original betas in:</p>
  
  <pre><code>&gt; f(betas)
  [1] 0
  </code></pre>
  
  <p>And relatively large with some random inputs:</p>
  
  <pre><code>&gt; set.seed(42)
  &gt; f(rnorm(5))
  [1] 8.849329
  </code></pre>
  
  <p>And minimized at the values I was able to calculate in Excel:</p>
  
  <pre><code>&gt; f(c(1.2,1.3,1.4,1.4,2.2))
  [1] 0.04
  </code></pre>