1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
|
# ifndef CPPAD_EXAMPLE_ABS_NORMAL_LP_BOX_HPP
# define CPPAD_EXAMPLE_ABS_NORMAL_LP_BOX_HPP
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
// SPDX-FileContributor: 2003-24 Bradley M. Bell
// ----------------------------------------------------------------------------
/*
{xrst_begin lp_box}
{xrst_spell
maxitr
rl
xout
}
abs_normal: Solve a Linear Program With Box Constraints
#######################################################
Syntax
******
| *ok* = ``lp_box`` (
| |tab| *level* , *A* , *b* , *c* , *d* , *maxitr* , *xout*
| )
Prototype
*********
{xrst_literal
// BEGIN PROTOTYPE
// END PROTOTYPE
}
Source
******
This following is a link to the source code for this example:
:ref:`lp_box.hpp-name` .
Problem
*******
We are given
:math:`A \in \B{R}^{m \times n}`,
:math:`b \in \B{R}^m`,
:math:`c \in \B{R}^n`,
:math:`d \in \B{R}^n`,
This routine solves the problem
.. math::
\begin{array}{rl}
\R{minimize} &
c^T x \; \R{w.r.t} \; x \in \B{R}^n
\\
\R{subject \; to} & A x + b \leq 0 \; \R{and} \; - d \leq x \leq d
\end{array}
Vector
******
The type *Vector* is a
simple vector with elements of type ``double`` .
level
*****
This value is less that or equal two.
If *level* == 0 ,
no tracing is printed.
If *level* >= 1 ,
a trace of the ``lp_box`` operations is printed.
If *level* >= 2 ,
the objective and primal variables :math:`x` are printed
at each :ref:`simplex_method-name` iteration.
If *level* == 3 ,
the simplex tableau is printed at each simplex iteration.
A
*
This is a :ref:`row-major<glossary@Row-major Representation>` representation
of the matrix :math:`A` in the problem.
b
*
This is the vector :math:`b` in the problem.
c
*
This is the vector :math:`c` in the problem.
d
*
This is the vector :math:`d` in the problem.
If :math:`d_j` is infinity, there is no limit for the size of
:math:`x_j`.
maxitr
******
This is the maximum number of newton iterations to try before giving up
on convergence.
xout
****
This argument has size is *n* and
the input value of its elements does no matter.
Upon return it is the primal variables
:math:`x` corresponding to the problem solution.
ok
**
If the return value *ok* is true, an optimal solution was found.
{xrst_toc_hidden
example/abs_normal/lp_box.cpp
example/abs_normal/lp_box.xrst
}
Example
*******
The file :ref:`lp_box.cpp-name` contains an example and test of
``lp_box`` .
{xrst_end lp_box}
-----------------------------------------------------------------------------
*/
# include "simplex_method.hpp"
// BEGIN C++
namespace CppAD { // BEGIN_CPPAD_NAMESPACE
// BEGIN PROTOTYPE
template <class Vector>
bool lp_box(
size_t level ,
const Vector& A ,
const Vector& b ,
const Vector& c ,
const Vector& d ,
size_t maxitr ,
Vector& xout )
// END PROTOTYPE
{ double inf = std::numeric_limits<double>::infinity();
//
size_t m = b.size();
size_t n = c.size();
//
CPPAD_ASSERT_KNOWN(
level <= 3, "lp_box: level is greater than 3");
CPPAD_ASSERT_KNOWN(
size_t(A.size()) == m * n, "lp_box: size of A is not m * n"
);
CPPAD_ASSERT_KNOWN(
size_t(d.size()) == n, "lp_box: size of d is not n"
);
if( level > 0 )
{ std::cout << "start lp_box\n";
CppAD::abs_print_mat("A", m, n, A);
CppAD::abs_print_mat("b", m, 1, b);
CppAD::abs_print_mat("c", n, 1, c);
CppAD::abs_print_mat("d", n, 1, d);
}
//
// count number of limits
size_t n_limit = 0;
for(size_t j = 0; j < n; j++)
{ if( d[j] < inf )
n_limit += 1;
}
//
// A_simplex and b_simplex define the extended constraints
Vector A_simplex((m + 2 * n_limit) * (2 * n) ), b_simplex(m + 2 * n_limit);
for(size_t i = 0; i < size_t(A_simplex.size()); i++)
A_simplex[i] = 0.0;
//
// put A * x + b <= 0 in A_simplex, b_simplex
for(size_t i = 0; i < m; i++)
{ b_simplex[i] = b[i];
for(size_t j = 0; j < n; j++)
{ // x_j^+ coefficient (positive component)
A_simplex[i * (2 * n) + 2 * j] = A[i * n + j];
// x_j^- coefficient (negative component)
A_simplex[i * (2 * n) + 2 * j + 1] = - A[i * n + j];
}
}
//
// put | x_j | <= d_j in A_simplex, b_simplex
size_t i_limit = 0;
for(size_t j = 0; j < n; j++) if( d[j] < inf )
{
// x_j^+ <= d_j constraint
b_simplex[ m + 2 * i_limit] = - d[j];
A_simplex[(m + 2 * i_limit) * (2 * n) + 2 * j] = 1.0;
//
// x_j^- <= d_j constraint
b_simplex[ m + 2 * i_limit + 1] = - d[j];
A_simplex[(m + 2 * i_limit + 1) * (2 * n) + 2 * j + 1] = 1.0;
//
++i_limit;
}
//
// c_simples
Vector c_simplex(2 * n);
for(size_t j = 0; j < n; j++)
{ // x_j+ component
c_simplex[2 * j] = c[j];
// x_j^- component
c_simplex[2 * j + 1] = - c[j];
}
size_t level_simplex = 0;
if( level >= 2 )
level_simplex = level - 1;
//
Vector x_simplex(2 * n);
bool ok = CppAD::simplex_method(
level_simplex, A_simplex, b_simplex, c_simplex, maxitr, x_simplex
);
for(size_t j = 0; j < n; j++)
xout[j] = x_simplex[2 * j] - x_simplex[2 * j + 1];
if( level > 0 )
{ CppAD::abs_print_mat("xout", n, 1, xout);
if( ok )
std::cout << "end lp_box: ok = true\n";
else
std::cout << "end lp_box: ok = false\n";
}
return ok;
}
} // END_CPPAD_NAMESPACE
// END C++
# endif
|
