We want to minimize the following objective function:

\[ f(x) = (x_1 - 1)^2 + (x_2 - 2)^2 \]

subject to the following constraints:

  1. Inequality constraint: \(g_1(x) = x_1^2 + x_2 - 1 \leq 0\)
  2. Equality constraint: \(h_1(x) = x_1 + x_2 - 1 = 0\)
# Install and load the necessary package
if (!require("nloptr")) {
  install.packages("nloptr")
}
## Loading required package: nloptr
# nloptr is an R package that serves as a bridge to NLopt, which is a free and open-source library designed for solving nonlinear 
# optimization problems. NLopt was initiated by Steven G. Johnson with the aim of creating a unified interface for accessing a wide range
# of available optimization routines scattered across the internet. In addition to these existing algorithms, NLopt also includes 
# original implementations of several other optimization techniques. Through the nloptr package, R users can conveniently access and 
# utilize the optimization capabilities offered by NLopt.
library(nloptr)

Define the objective function

# Define the objective function
objective_function <- function(x) {
  return((x[1] - 1)^2 + (x[2] - 2)^2)
}

Define the gradient of the objective function

# Define the gradient of the objective function
gradient_objective_function <- function(x) {
  grad <- c(2 * (x[1] - 1), 2 * (x[2] - 2))
  return(grad)
}

Define the inequality constraint

# Define the inequality constraint
inequality_constraint <- function(x) {
  return(x[1]^2 + x[2] - 1)
}

Define the gradient of the inequality constraint

# Define the gradient of the inequality constraint
gradient_inequality_constraint <- function(x) {
  grad <- c(2 * x[1], 1)
  return(grad)
}

Define the equality constraint

# Define the equality constraint
equality_constraint <- function(x) {
  return(x[1] + x[2] - 1)
}

Define the gradient of the equality constraint

# Define the gradient of the equality constraint
gradient_equality_constraint <- function(x) {
  grad <- c(1, 1)
  return(grad)
}

Initial guess for the variables

# Initial guess for the variables
x0 <- c(0.5, 0.5)

Definition of the optimization problem

# Define the optimization problem using NLOPT_LD_AUGLAG with NLOPT_LD_LBFGS as a sub-algorithm
result <- nloptr(
  x0 = x0,
  eval_f = objective_function,
  eval_grad_f = gradient_objective_function,
  lb = c(-Inf, -Inf),
  ub = c(Inf, Inf),
  eval_g_ineq = function(x) { return(c(inequality_constraint(x))) },
  eval_jac_g_ineq = function(x) { return(matrix(gradient_inequality_constraint(x), nrow = 1)) },
  eval_g_eq = function(x) { return(c(equality_constraint(x))) },
  eval_jac_g_eq = function(x) { return(matrix(gradient_equality_constraint(x), nrow = 1)) },
  opts = list(
    algorithm = "NLOPT_LD_AUGLAG",
    local_opts = list(
      algorithm = "NLOPT_LD_LBFGS",
      xtol_rel = 1.0e-8
    ),
    xtol_rel = 1.0e-8
  )
)

Solve the optimization problem and print the result

result <- nloptr(
  x0 = x0,
  eval_f = objective_function,
  eval_grad_f = gradient_objective_function,
  lb = c(-Inf, -Inf),
  ub = c(Inf, Inf),
  eval_g_ineq = function(x) { return(c(inequality_constraint(x))) },
  eval_jac_g_ineq = function(x) { return(matrix(gradient_inequality_constraint(x), nrow = 1)) },
  eval_g_eq = function(x) { return(c(equality_constraint(x))) },
  eval_jac_g_eq = function(x) { return(matrix(gradient_equality_constraint(x), nrow = 1)) },
  opts = list(
    algorithm = "NLOPT_LD_AUGLAG",
    local_opts = list(
      algorithm = "NLOPT_LD_LBFGS",
      xtol_rel = 1.0e-8
    ),
    xtol_rel = 1.0e-8,
    maxeval = 1000 # Increase the maximum number of evaluations
  )
)


# Print the result
print(result)
## 
## Call:
## 
## nloptr(x0 = x0, eval_f = objective_function, eval_grad_f = gradient_objective_function, 
##     lb = c(-Inf, -Inf), ub = c(Inf, Inf), eval_g_ineq = function(x) {
##         return(c(inequality_constraint(x)))
##     }, eval_jac_g_ineq = function(x) {
##         return(matrix(gradient_inequality_constraint(x), nrow = 1))
##     }, eval_g_eq = function(x) {        return(c(equality_constraint(x)))
##     }, eval_jac_g_eq = function(x) {
##         return(matrix(gradient_equality_constraint(x), nrow = 1))
##     }, opts = list(algorithm = "NLOPT_LD_AUGLAG", local_opts = list(algorithm = "NLOPT_LD_LBFGS", 
##         xtol_rel = 1e-08), xtol_rel = 1e-08, maxeval = 1000))
## 
## 
## Minimization using NLopt version 2.7.1 
## 
## NLopt solver status: 4 ( NLOPT_XTOL_REACHED: Optimization stopped because 
## xtol_rel or xtol_abs (above) was reached. )
## 
## Number of Iterations....: 188 
## Termination conditions:  xtol_rel: 1e-08 maxeval: 1000 
## Number of inequality constraints:  1 
## Number of equality constraints:    1 
## Optimal value of objective function:  1.99999999949347 
## Optimal value of controls: 3.694613e-08 1
library(ggplot2)
# Plot the objective function
x1 <- seq(-1, 3, length.out = 100)
x2 <- seq(0, 4, length.out = 100)
f <- outer(x1, x2, function(x1, x2) (x1 - 1)^2 + (x2 - 2)^2)

df <- expand.grid(x1 = x1, x2 = x2)
df$z <- as.vector(f)

ggplot(df, aes(x = x1, y = x2, z = z)) +
  geom_contour_filled() +
  labs(title = "Contour Plot of the Objective Function",
       x = expression(x[1]),
       y = expression(x[2])) +
  theme_minimal()

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
# Plot the objective function in 3D
x1 <- seq(-1, 3, length.out = 100)
x2 <- seq(0, 4, length.out = 100)
f <- outer(x1, x2, function(x1, x2) (x1 - 1)^2 + (x2 - 2)^2)

plot_ly(x = ~x1, y = ~x2, z = ~f) %>%
  add_surface() %>%
  layout(
    title = "3D Plot of the Objective Function",
    scene = list(
      xaxis = list(title = "x1"),
      yaxis = list(title = "x2"),
      zaxis = list(title = "f(x)")
    )
  )