# source:stable/3.9/Ipopt/contrib/RInterface/tests/hs071_nlp.R@1831

Last change on this file since 1831 was 1831, checked in by andreasw, 3 years ago

merged with trunk rev 1830

File size: 3.6 KB
Line
2# This code is published under the Common Public License.
3#
4# File:   hs071_nlp.R
5# Author: Jelmer Ypma
6# Date:   18 April 2010
7#
8# Example problem, number 71 from the Hock-Schittkowsky test suite
9#
10# \min_{x} x1*x4*(x1 + x2 + x3) + x3
11# s.t.
12#    x1*x2*x3*x4 >= 25
13#    x1^2 + x2^2 + x3^2 + x4^2 = 40
14#    1 <= x1,x2,x3,x4 <= 5
15#
16# x0 = (1,5,5,1)
17#
18# optimal solution = (1.00000000, 4.74299963, 3.82114998, 1.37940829)
19#
20# Adapted from the Ipopt C++ interface example.
21
22library('ipoptr')
23
24#
25# f(x) = x1*x4*(x1 + x2 + x3) + x3
26#
27eval_f <- function( x ) {
28    return( x[1]*x[4]*(x[1] + x[2] + x[3]) + x[3] )
29}
30
31eval_grad_f <- function( x ) {
32    return( c( x[1] * x[4] + x[4] * (x[1] + x[2] + x[3]),
33               x[1] * x[4],
34               x[1] * x[4] + 1.0,
35               x[1] * (x[1] + x[2] + x[3]) ) )
36}
37
38# constraint functions
39eval_g <- function( x ) {
40    return( c( x[1] * x[2] * x[3] * x[4],
41               x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2 ) )
42}
43
44# The Jacobian for this problem is dense
45eval_jac_g_structure <- list( c(1,2,3,4), c(1,2,3,4) )
46
47eval_jac_g <- function( x ) {
48    return( c ( x[2]*x[3]*x[4],
49                x[1]*x[3]*x[4],
50                x[1]*x[2]*x[4],
51                x[1]*x[2]*x[3],
52                2.0*x[1],
53                2.0*x[2],
54                2.0*x[3],
55                2.0*x[4] ) )
56}
57
58# The Hessian for this problem is actually dense,
59# This is a symmetric matrix, fill the lower left triangle only.
60eval_h_structure <- list( c(1), c(1,2), c(1,2,3), c(1,2,3,4) )
61
62eval_h <- function( x, obj_factor, hessian_lambda ) {
63
64    values <- numeric(10)
65    values[1] = obj_factor * (2*x[4]) # 1,1
66
67    values[2] = obj_factor * (x[4])   # 2,1
68    values[3] = 0                     # 2,2
69
70    values[4] = obj_factor * (x[4])   # 3,1
71    values[5] = 0                     # 4,2
72    values[6] = 0                     # 3,3
73
74    values[7] = obj_factor * (2*x[1] + x[2] + x[3]) # 4,1
75    values[8] = obj_factor * (x[1])                 # 4,2
76    values[9] = obj_factor * (x[1])                 # 4,3
77    values[10] = 0                                  # 4,4
78
79
80    # add the portion for the first constraint
81    values[2] = values[2] + hessian_lambda[1] * (x[3] * x[4]) # 2,1
82
83    values[4] = values[4] + hessian_lambda[1] * (x[2] * x[4]) # 3,1
84    values[5] = values[5] + hessian_lambda[1] * (x[1] * x[4]) # 3,2
85
86    values[7] = values[7] + hessian_lambda[1] * (x[2] * x[3]) # 4,1
87    values[8] = values[8] + hessian_lambda[1] * (x[1] * x[3]) # 4,2
88    values[9] = values[9] + hessian_lambda[1] * (x[1] * x[2]) # 4,3
89
90    # add the portion for the second constraint
91    values[1] = values[1] + hessian_lambda[2] * 2 # 1,1
92    values[3] = values[3] + hessian_lambda[2] * 2 # 2,2
93    values[6] = values[6] + hessian_lambda[2] * 2 # 3,3
94    values[10] = values[10] + hessian_lambda[2] * 2 # 4,4
95
96    return ( values )
97}
98
99# initial values
100x0 <- c( 1, 5, 5, 1 )
101
102# lower and upper bounds of control
103lb <- c( 1, 1, 1, 1 )
104ub <- c( 5, 5, 5, 5 )
105
106# lower and upper bounds of constraints
107constraint_lb <- c(  25, 40 )
108constraint_ub <- c( Inf, 40 )
109
110
111opts <- list("print_level"=0,
112             "file_print_level"=12,
113             "output_file"="hs071_nlp.out")
114
115print( ipoptr( x0=x0,
116               eval_f=eval_f,