#はじめに
Rで最適化問題を解くことができるらしいので，学部レベル(初級)のミクロ経済学の序盤で習う，「予算制約の中で効用（満足感）を最大化する購買量の探索」という効用最大化問題（最適消費計画問題）を例にとり，Rで簡単な制約付き最適化問題を解いたりイロイロしてみます。

◇2018年1月20日追記

ここに載せたことをshinyで作動させてみました。
シンプルな効用最大化問題計算機という感じです。
https://nigimitama.shinyapps.io/MicroEconomics_01_UtilityMaximizationProblem/

中身のコードはこのページの下部（おまけ2）にあります。よろしければ御覧ください。

#問題
2財モデルのこんな問題があるとします。

第1財の消費量を$x_1$，第2財の消費量を$x_2$と表す。ある家計の効用関数が$U(x_1,x_2)=x_1^{\frac{2}{5}}x_2^{\frac{3}{5}}$で表されるとする。予算が$100$，第1財の価格が$4$，第2財の価格が$6$のとき，効用を最大化する最適消費量はそれぞれいくらか？

$$\max_{x_1,x_2} \hspace{1em} x_1^{\frac{2}{5}}x_2^{\frac{3}{5}}$$
$$s.t. \hspace{1em} 100 = 4x_1+6x_2$$

ラグランジュの未定乗数法で解ける問題ですが，Rで解いてみます。

#Rで解く
{Rsolnp}パッケージを使います。

library(Rsolnp)
ObjFunc = function(x) return( - x[1]^(2/5) * x[2]^(3/5) ) #目的関数
#solnp()は最小化をするようになっているので，目的関数にマイナスを掛ける

ConstFunc = function(x) return( x[1]* 4 + x[2] * 6 ) #制約式の右辺
eq.value <- c(100) #制約式の左辺

x0 <- c(1,1) #決定変数を初期化

solution <- solnp(x0, fun = ObjFunc, eqfun = ConstFunc, eqB = eq.value )

# 最適な｢x_1｣と「x_2｣の値
solution$pars

[1] 10 10

最適な消費量$x_1^* , x_2^* $ は $x_1^* = 10, \ x_2^* = 10$であることが示された

#（おまけ1）図示
####■効用関数の図示
モノ（$x_1,x_2$）をたくさん買うほど効用$U$は上がる

# 3次元プロット:persp()
x_1 <- 1:50
x_2 <- 1:50
u <- function(x_1,x_2) {x_1^(2/5) * x_2^(3/5)} #効用関数を定義
U <- outer(x_1, x_2, u) #outer()はx_1,x_2に対応したf(x_1,x_2)の値を行列で返す
persp(x_1, x_2, U,
      theta = 30, # 横回転の角度
      phi = 30, # 縦回転の角度
      ticktype = "simple", # 線の種類
      lwd = 0.5, # 線の太さ
      col = F,
      border = 8)

####■等高線（無差別曲線）と予算制約線と最適消費点
効用関数の曲面を真上から見て等高線にした図。経済学ではこの等高線を無差別曲線という。
予算制約線（青い直線）と最適消費点（藍色の点）も記しておいた。

無差別曲線と予算制約線が接する点が最適消費点となる。

# 等高線:contour()
contour(x_1, x_2, U, method = "edge", labcex = 1,lwd = 2)
abline(a = 100/6, b = -4/6, lwd = 2, col = "blue") #予算制約線
points(x = 10, y = 10, lwd = 3, col = "darkblue", pch = 16) #最適消費点

#（おまけ2）Shinyで効用最大化問題計算機を作る
以上の制約付き最適化問題を解く{Rsolnp}の処理や関数の図示などをshinyで行わせてみる。

ui.R}

# 日本語の文字列が本文中にあるとdeployできないみたい
# server側のRenderUI内の日本語は問題ないのでコレを利用

library(shinydashboard)
library(shiny)
library(Rsolnp)

header <- dashboardHeader(
  title = "Utility maximization problem"
)

body <- dashboardBody(
    column(width = 3,
           #予算制約
           box(width = NULL,
               uiOutput('budgettext'),
               numericInput("budget", 
                            label = "$$B$$", 
                            value = 100),
               numericInput("p_1", 
                            label = "$$p_1$$", 
                            value = 2),
               numericInput("p_2", 
                            label = "$$p_2$$", 
                            value = 2),
               uiOutput("budgetform")
              
           ),
           
           box(width = NULL, 
               #効用関数モジュール
               uiOutput('utilityfunctiontext'),
               withMathJax(),   #TeX表記有効化
               numericInput("a", 
                            label = "$$a$$", 
                            value = 1),
               numericInput("b", 
                            label = "$$b$$", 
                            value = 1),
               numericInput("c", 
                            label = "$$c$$", 
                            value = 1),
               numericInput("d", 
                            label = "$$d$$", 
                            value = 2),
               numericInput("e", 
                            label = "$$e$$", 
                            value = 1),
               numericInput("f", 
                            label = "$$f$$", 
                            value = 2),
               uiOutput("utilityfunction")
               ),
           box(width = NULL,
               #p("変数の範囲（描画範囲）"),
               numericInput("range", 
                            label =  "range of graph", 
                            value = 50)
               )
    ),
    fluidRow(
      column(width = 7,
             #設定された効用最大化問題
             box(width = NULL,
                 uiOutput("Utilitymaximizationproblem")
             ),
             #効用最大化問題の解
             box(width = NULL, status = "warning",
                 uiOutput("maxutility")
                 ),
             #無差別曲線
             box(width = NULL, solidHeader = TRUE,
                 uiOutput('indifferencecurve'),
                 plotOutput("Plot2D", width = 300, height = 300)
                 ),
             box(width = NULL, solidHeader = TRUE,
                 uiOutput('Plot3Dtext'),
                 plotOutput("Plot3D",height = 300)
                 )
             )
      )
)

dashboardPage(
  header,
  dashboardSidebar(disable = TRUE),
  body
)

server.R}

library(shinydashboard)
library(shiny)
library(Rsolnp)

# サーバロジックの定義
shinyServer(function(input, output) {
  
  # 予算制約式box見出し
  output$budgettext <- renderUI({
    withMathJax(helpText("予算制約と価格の設定$$p_1 x_1 + p_2 x_2 \\leq B$$"))
  })
  #指定された予算をTeXでUIに表示
  output$budgetform <- renderUI({
    withMathJax(helpText(paste0('指定された予算制約式:$$',
                                input$p_1,'x_1 + ', input$p_2,'x_2',
                                '\\leq',input$budget,'$$')))
  })
  
  
  # 効用関数box見出し
  output$utilityfunctiontext <- renderUI({
    withMathJax(helpText('効用関数のパラメータを設定$$U(x_1,x_2)=a x_1^ \\frac{c}{d} b x_2^ \\frac{e}{f}$$'))
  })
  #指定された効用関数をTeXでUIに表示
  output$utilityfunction <- renderUI({
    withMathJax(helpText(paste0('指定された効用関数:$$U(x_1,x_2)=',
                              input$a,'x_1','^\\frac{',input$c,'}','{',input$d,'}', input$b,'x_2','^\\frac{',input$e,'}','{',input$f,'}','$$')))
  })
  
  
  #設定された効用最大化問題
  output$Utilitymaximizationproblem <- renderUI({
    withMathJax(helpText(paste0(
      '設定された効用最大化問題:',

      '$$ \\max_{x_1,x_2} \\hspace{1em}  U(x_1,x_2)=',input$a,'x_1','^\\frac{',input$c,'}','{',input$d,'}', input$b,'x_2','^\\frac{',input$e,'}','{',input$f,'}$$',
      
      '$$s.t. \\hspace{1em} ', input$p_1,'x_1 + ', input$p_2,'x_2 \\leq',input$budget,'$$'
      )))
  })
  
  
  #最適消費計画の解
  output$maxutility <- renderUI({
    
    ObjFunc = function(x) return( - input$a * x[1]^(input$c/input$d) * input$b * x[2]^(input$e/input$f) ) #目的関数
    #solnp()は最小化をするようになっているので，目的関数にマイナスを掛ける
    ConstFunc = function(x) return( x[1]* input$p_1 + x[2] * input$p_2 ) #制約式の右辺
    eq.value <- c(input$budget) #制約式の左辺
    x0 <- c(1,1) #決定変数を初期化
    
    solution <- solnp(x0, fun = ObjFunc, eqfun = ConstFunc, eqB = eq.value )
    
    #return
    withMathJax(helpText(paste0("効用を最大化する消費量:",
                                '$$x_{1}^*=',solution$pars[1],', \\hspace{1em}  x_{2}^*=',solution$pars[2],'$$')))
  })
  
  
  #無差別曲線と予算制約線と最適消費点
  output$indifferencecurve <- renderUI({
    withMathJax(helpText('無差別曲線と予算制約線と最適消費点'))
  })
  output$Plot2D <- renderPlot({
    #以下のreactiveな要素をrenderPlotに入れれば使えるっぽい
    #効用関数をfunctionの形で定義
    u <- function(x_1,x_2) {input$a * x_1 ^ (input$c/input$d) * input$b * x_2 ^ (input$e/input$f)} #効用関数を定義
    
    #xの範囲指定
    x_1 <- 1:input$range
    x_2 <- 1:input$range
    U <- outer(x_1, x_2, u) #outer()はx_1,x_2に対応したf(x_1,x_2)の値を行列で返す
    
    #plot
    par(mar=c(3,3,1,1))
    contour(x_1, x_2, U, method = "edge", labcex = 1,lwd = 2)
    abline(a = input$budget/input$p_2, b = -input$p_1/input$p_2, lwd = 2, col = "blue") #予算制約線 x_2 = ...
    
    #最適消費点
    ObjFunc = function(x) return( - input$a * x[1]^(input$c/input$d) * input$b * x[2]^(input$e/input$f) ) #目的関数
    #solnp()は最小化をするようになっているので，目的関数にマイナスを掛ける
    ConstFunc = function(x) return( x[1]* input$p_1 + x[2] * input$p_2 ) #制約式の右辺
    eq.value <- c(input$budget) #制約式の左辺
    x0 <- c(1,1) #決定変数を初期化
    solution <- solnp(x0, fun = ObjFunc, eqfun = ConstFunc, eqB = eq.value )
    points(x = solution$pars[1], y = solution$pars[2], lwd = 3, col = "darkblue", pch = 16) #最適消費点
    
  })  
  
  #3D plot of utility function 効用関数の3次元図
  output$Plot3Dtext <- renderUI({
    withMathJax(helpText('効用関数の3次元図'))
  })
  output$Plot3D <- renderPlot({
    
    #以下のreactiveな要素をrenderPlotに入れれば使えるっぽい
    #効用関数をfunctionの形で定義
    u <- function(x_1,x_2) {input$a * x_1 ^ (input$c/input$d) * input$b * x_2 ^ (input$e/input$f)} #効用関数を定義
    
    #xの範囲指定
    x_1 <- 1:input$range
    x_2 <- 1:input$range
    U <- outer(x_1, x_2, u) #outer()はx_1,x_2に対応したf(x_1,x_2)の値を行列で返す    
    persp(x_1, x_2, U,
          theta = 30, # 横回転の角度。ここをいじって表示のアングルを変える。
          phi = 20, # 縦回転の角度。ここをいじって表示のアングルを変える。
          ticktype = "simple", # 線の種類
          lwd = 0.5, # 線が太いと見づらかったりするので0.5にしておいた。
          col = F, #塗りつぶしなし
          border = 8 #枠線灰色
          )
  })  
})