1% Copyright (C) 2009 International Business Machines
2% All Rights Reserved.
3% This code is published under the Eclipse Public License.
4%
5% $Id: hs071_c.c 699 2006-04-05 21:05:18Z andreasw $
6%
7% Author:  Andreas Waechter               IBM    2009-04-02
8%
9% This file is part of the Ipopt tutorial.  It is a version with
10% mistakes for the matlab implemention of the coding exercise problem
11% (in AMPL formulation):
12%
13% param n := 4;
14%
15% var x {1..n} <= 0, >= -1.5, := -0.5;
16%
17% minimize obj:
18%   sum{i in 1..n} (x[i]-1)^2;
19%
20% subject to constr {i in 2..n-1}:
21%   (x[i]^2+1.5*x[i]-i/n)*cos(x[i+1]) - x[i-1] = 0;
22%
23% The constant term "i/n" in the constraint is supposed to be input data
24
25function [x, info] = TutorialMatlab
26
27  % Size of the problem
28  n = 5;
29
30  % Problem data
31  a = ((1:n)/n)';
32
33  % Starting point
34  x0 = -0.5*ones(n,1);
35
36  % Lower and upper bounds for the variables
37  options.lb = -1.5*ones(n,1);
38  options.ub = zeros(n,1);
39
40  % Constraint bounds
41  options.cl = zeros(n-2,1);
42  options.cu = zeros(n-2,1);
43
44  % Set the Ipopt options
45  % options.ipopt.mu_strategy = 'adaptive';
46  options.ipopt.derivative_test = 'second-order';
47  options.ipopt.max_iter = 10;
48
49  % Set the callback functions
50  funcs.objective         = @eval_f;
51  funcs.constraints       = @eval_g;
52  funcs.gradient          = @eval_grad_f;
53  funcs.jacobian          = @eval_jac_g;
54  funcs.jacobianstructure = @eval_jac_g_struct;
55  funcs.hessian           = @eval_hess;
56  funcs.hessianstructure  = @eval_hess_struct;
57
58  [x info] = ipopt(x0, funcs, options);
59
60
61  % End of main function
62  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63
64  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
65  % Evaluate value of objective function
66  function f = eval_f(x)
67
68    tmp = (x - 1).^2;
69    f = sum(tmp);
70
71  end
72
73  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74  % Evaluate gradient of objective function
75  function df = eval_grad_f(x)
76
77    df = (x - 1);
78
79  end
80
81  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82  % Evaluate value of constraint bodies
83  function g = eval_g(x)
84
85    g = (x(2:n-1).^2 + 1.5*x(1:n-2) - a(2:n-1)).*cos(x(3:n)) - x(1:n-2);
86
87  end
88
89  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
90  % Return constraint Jacobian strcture
91  function A = eval_jac_g_struct
92
93    % tri-diagonal structure
94    Diags = ones(n, 3);
95    A = spdiags(Diags, 0:2, n-2, n);
96
97  end
98
99  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100  % Evaluate constraint Jacobian
101  function A = eval_jac_g(x)
102
103    Diags = -1*ones(n-2,3);
104    Diags(:,2) = (2*x(2:n-1)+1.5).*cos(x(3:n));
105    Diags(:,3) = -(x(2:n-1).^2 + 1.5*x(2:n-1) - a(2:n-1)).*sin(x(3:n));
106    A = spdiags(Diags, 0:2, n-2, n);
107
108  end
109
110  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
111  % Return Hessian of Lagrangian function structure
112  function H = eval_hess_struct
113
114    Diags = ones(n, 2);
115    H = spdiags(Diags, -1:0, n, n);
116
117  end
118
119  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
120  % Evaluate Hessian of Lagrangian function
121  function H = eval_hess(x, sigma, lambda)
122
123    Diags = zeros(n,2);
124
125    % part from the objective function
126    Diags(:,2) = sigma * 2;
127
128    % (x_i , x_i) part
129    Diags(2:n-1,2) = Diags(2:n-1,2) + lambda.*cos(x(3:n));
130
131    % (x_{i+1}, x_{i+1}) part
132    Diags(3:n,2) = Diags(3:n,2) ...
133      - lambda.*(x(2:n-1).^2 + 1.5*x(1:n-2) - a(2:n-1)).*cos(x(3:n));
134
135    % (x_i, x_{i+1}) part
136    Diags(2:n-1,1) = -lambda.*sin(x(3:n)).*(2*x(2:n-1)+1.5);
137
138    H = spdiags(Diags, -1:0, n, n);
139
140  end
141
142end
143