* Bug #14582: gettext or _ were applied to broken literal strings
[scilab.git] / scilab / modules / optimization / help / en_US / karmarkar.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <!--
3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4  * Copyright (C) 2008 - INRIA
5  * Copyright (C) 2010 - 2011 - DIGITEO - Michael Baudin
6  * 
7  * Copyright (C) 2012 - 2016 - Scilab Enterprises
8  *
9  * This file is hereby licensed under the terms of the GNU GPL v2.0,
10  * pursuant to article 5.3.4 of the CeCILL v.2.1.
11  * This file was originally licensed under the terms of the CeCILL v2.1,
12  * and continues to be available under such terms.
13  * For more information, see the COPYING file which you should have received
14  * along with this program.
15  *
16  -->
17 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="karmarkar" xml:lang="en">
18     <refnamediv>
19         <refname>karmarkar</refname>
20         <refpurpose>Solves a linear optimization problem.</refpurpose>
21     </refnamediv>
22     <refsynopsisdiv>
23         <title>Syntax</title>
24         <synopsis>
25             xopt=karmarkar(Aeq,beq,c)
26             xopt=karmarkar(Aeq,beq,c,x0)
27             xopt=karmarkar(Aeq,beq,c,x0,rtolf)
28             xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam)
29             xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter)
30             xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun)
31             xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun,A,b)
32             xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun,A,b,lb)
33             xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun,A,b,lb,ub)
34             [xopt,fopt] = karmarkar(...)
35             [xopt,fopt,exitflag] = karmarkar(...)
36             [xopt,fopt,exitflag,iter] = karmarkar(...)
37             [xopt,fopt,exitflag,iter,yopt] = karmarkar(...)
38         </synopsis>
39     </refsynopsisdiv>
40     <refsection>
41         <title>Arguments</title>
42         <variablelist>
43             <varlistentry>
44                 <term>Aeq</term>
45                 <listitem>
46                     <para>
47                         a n-by-p matrix of doubles, where <literal>n</literal> is the number
48                         of linear equality constraints
49                         and <literal>p</literal> is the number of unknowns,
50                         the matrix in the linear equality constraints.
51                     </para>
52                 </listitem>
53             </varlistentry>
54             <varlistentry>
55                 <term>beq</term>
56                 <listitem>
57                     <para>
58                         a n-by-1 matrix of doubles,
59                         the right hand side of the linear equality constraint.
60                     </para>
61                 </listitem>
62             </varlistentry>
63             <varlistentry>
64                 <term>c</term>
65                 <listitem>
66                     <para>
67                         a p-by-1 matrix of doubles,
68                         the linear part of the objective function.
69                     </para>
70                 </listitem>
71             </varlistentry>
72             <varlistentry>
73                 <term>x0</term>
74                 <listitem>
75                     <para>
76                         a p-by-1 matrix of doubles, the initial guess (default <literal>x0=[]</literal>).
77                         If <literal>x0==[]</literal>, the karmarkar function automatically computes a feasible initial
78                         guess.
79                     </para>
80                 </listitem>
81             </varlistentry>
82             <varlistentry>
83                 <term>rtolf</term>
84                 <listitem>
85                     <para>
86                         a 1-by-1 matrix of doubles,
87                         a relative tolerance on <literal>f(x)=c'*x</literal> (default <literal>rtolf=1.d-5</literal>).
88                     </para>
89                 </listitem>
90             </varlistentry>
91             <varlistentry>
92                 <term>gam</term>
93                 <listitem>
94                     <para>
95                         a 1-by-1 matrix of doubles, the scaling factor (default <literal>gam=0.5</literal>).
96                         The scaling factor must satisfy <literal>0&lt;gam&lt;1</literal>.
97                     </para>
98                 </listitem>
99             </varlistentry>
100             <varlistentry>
101                 <term>maxiter</term>
102                 <listitem>
103                     <para>
104                         a 1-by-1 matrix of floating point integers, the maximum number of iterations (default <literal>maxiter=200</literal>).
105                         The maximum number of iterations must be greater than 1.
106                     </para>
107                 </listitem>
108             </varlistentry>
109             <varlistentry>
110                 <term>outfun</term>
111                 <listitem>
112                     <para>
113                         a function or a list, the output function. See below for details.
114                     </para>
115                 </listitem>
116             </varlistentry>
117             <varlistentry>
118                 <term>A</term>
119                 <listitem>
120                     <para>
121                         a ni-by-p matrix of doubles, the matrix of linear inequality constraints.
122                     </para>
123                 </listitem>
124             </varlistentry>
125             <varlistentry>
126                 <term>b</term>
127                 <listitem>
128                     <para>
129                         a ni-by-1 matrix of doubles, the right-hand side of linear inequality constraints.
130                     </para>
131                 </listitem>
132             </varlistentry>
133             <varlistentry>
134                 <term>lb</term>
135                 <listitem>
136                     <para>
137                         a p-by-1 matrix of doubles, the lower bounds.
138                     </para>
139                 </listitem>
140             </varlistentry>
141             <varlistentry>
142                 <term>ub</term>
143                 <listitem>
144                     <para>
145                         a p-by-1 matrix of doubles, the upper bounds.
146                     </para>
147                 </listitem>
148             </varlistentry>
149             <varlistentry>
150                 <term>xopt</term>
151                 <listitem>
152                     <para>a p-by-1 matrix of doubles, the solution.</para>
153                 </listitem>
154             </varlistentry>
155             <varlistentry>
156                 <term>fopt</term>
157                 <listitem>
158                     <para>
159                         a 1-by-1 matrix of doubles, the objective function
160                         value at optimum, i.e. <literal>fopt=c'*xopt</literal>.
161                     </para>
162                 </listitem>
163             </varlistentry>
164             <varlistentry>
165                 <term>exitflag</term>
166                 <listitem>
167                     <para>
168                         a 1-by-1 matrix of floating point integers, the status of execution.
169                         See below for details.
170                     </para>
171                 </listitem>
172             </varlistentry>
173             <varlistentry>
174                 <term>iter</term>
175                 <listitem>
176                     <para>
177                         a 1-by-1 matrix of floating point integers, the number of iterations.
178                     </para>
179                 </listitem>
180             </varlistentry>
181             <varlistentry>
182                 <term>yopt</term>
183                 <listitem>
184                     <para>
185                         a <literal>struct</literal> containing the dual solution.
186                         The structure yopt has four fields : ineqlin, eqlin, upper, lower.
187                         The field <literal>yopt.ineqlin</literal> is the Lagrange multiplier for the inequality constraints.
188                         The field <literal>yopt.eqlin</literal> is the Lagrange multiplier for the equality constraints.
189                         The field <literal>yopt.upper</literal> is the Lagrange multiplier for the upper bounds.
190                         The field <literal>yopt.lower</literal> is the Lagrange multiplier for the lower bounds.
191                     </para>
192                 </listitem>
193             </varlistentry>
194         </variablelist>
195     </refsection>
196     <refsection>
197         <title>Description</title>
198         <para>
199             Computes the solution of linear programming problems.
200         </para>
201         <para>
202             This function has the two following modes.
203         </para>
204         <itemizedlist>
205             <listitem>
206                 <para>
207                     If no inequality constraints and no bound is given
208                     (i.e. if <literal>(A==[] &amp; lb==[] &amp; ub==[])</literal>),
209                     the function ensures that the variable is nonnegative.
210                 </para>
211             </listitem>
212             <listitem>
213                 <para>
214                     If any inequality constraints or any bound is given
215                     (i.e. if <literal>(A&lt;&gt;[] or lb&lt;&gt;[] or ub&lt;&gt;[])</literal>),
216                     the function takes into account for this inequality or bound
217                     (and does not ensure that the variable is nonnegative).
218                 </para>
219             </listitem>
220         </itemizedlist>
221         <para>
222             If no inequality constraints and no bound is given
223             (i.e. if <literal>(A==[] &amp; lb==[] &amp; ub==[])</literal>),
224             solves the linear optimization problem:
225         </para>
226         <para>
227             <latex>
228                 \begin{eqnarray}
229                 \begin{array}{l}
230                 \textrm{minimize } c^T \cdot x\\
231                 A_{eq} x = b_{eq}\\
232                 x \geq 0
233                 \end{array}
234                 \end{eqnarray}
235             </latex>
236         </para>
237         <para>
238             If any inequality constraints or any bound is given
239             (i.e. if <literal>(A&lt;&gt;[] | lb&lt;&gt;[] | ub&lt;&gt;[])</literal>),
240             solves the linear optimization problem:
241         </para>
242         <para>
243             <latex>
244                 \begin{eqnarray}
245                 \begin{array}{l}
246                 \textrm{minimize } c^T \cdot x\\
247                 A_{eq} x = b_{eq}\\
248                 A x \leq b\\
249                 \ell_b \leq x \leq u_b
250                 \end{array}
251                 \end{eqnarray}
252             </latex>
253         </para>
254         <para>
255             Any optional parameter equal to the empty matrix <literal>[]</literal> is replaced by
256             its default value.
257         </para>
258         <para>
259             The <literal>exitflag</literal> parameter allows to know why the algorithm
260             terminated.
261         </para>
262         <itemizedlist>
263             <listitem>
264                 <para>
265                     <literal>exitflag = 1</literal> if algorithm converged.
266                 </para>
267             </listitem>
268             <listitem>
269                 <para>
270                     <literal>exitflag = 0</literal> if maximum number of iterations was reached.
271                 </para>
272             </listitem>
273             <listitem>
274                 <para>
275                     <literal>exitflag = -1</literal> if no feasible point was found
276                 </para>
277             </listitem>
278             <listitem>
279                 <para>
280                     <literal>exitflag = -2</literal> if problem is unbounded.
281                 </para>
282             </listitem>
283             <listitem>
284                 <para>
285                     <literal>exitflag = -3</literal> if search direction became zero.
286                 </para>
287             </listitem>
288             <listitem>
289                 <para>
290                     <literal>exitflag = -4</literal> if algorithm stopped on user's request.
291                 </para>
292             </listitem>
293         </itemizedlist>
294         <para>
295             The output function <literal>outfun</literal> must have header
296         </para>
297         <programlisting role="example"><![CDATA[ 
298 stop = outfun ( xopt , optimValues , state )
299 ]]></programlisting>
300         <para>
301             where
302             <literal>xopt</literal> is a p-by-1 matrix of double representing
303             the current point,
304             <literal>optimValues</literal> is a <literal>struct</literal>,
305             <literal>state</literal> is a 1-by-1 matrix of strings.
306             Here, <literal>stop</literal> is a 1-by-1 matrix of booleans,
307             which is <literal>%t</literal> if the algorithm must stop.
308         </para>
309         <para>
310             <literal>optimValues</literal> has the following fields:
311         </para>
312         <itemizedlist>
313             <listitem>
314                 <para>
315                     <literal>optimValues.funccount</literal>: a 1-by-1 matrix of floating point integers,
316                     the number of function evaluations
317                 </para>
318             </listitem>
319             <listitem>
320                 <para>
321                     <literal>optimValues.fval</literal>: a 1-by-1 matrix of doubles, the best function value
322                 </para>
323             </listitem>
324             <listitem>
325                 <para>
326                     <literal>optimValues.iteration</literal>: a 1-by-1 matrix of floating point integers,
327                     the current iteration number
328                 </para>
329             </listitem>
330             <listitem>
331                 <para>
332                     <literal>optimValues.procedure</literal>: a 1-by-1 matrix of strings, the type of step performed.
333                 </para>
334             </listitem>
335             <listitem>
336                 <para>
337                     <literal>optimValues.dualgap</literal>: a 1-by-1 matrix of doubles, the duality gap, i.e.
338                     <literal>abs(yopt'*beq - fopt)</literal>.
339                     At optimum, the duality gap is zero.
340                 </para>
341             </listitem>
342         </itemizedlist>
343         <para>
344             The <literal>optimValues.procedure</literal> field can have the following values.
345         </para>
346         <itemizedlist>
347             <listitem>
348                 <para>
349                     If <literal>optimValues.procedure="x0"</literal>, then the algorithm is computing the initial
350                     feasible guess <literal>x0</literal> (phase 1).
351                 </para>
352             </listitem>
353             <listitem>
354                 <para>
355                     If <literal>optimValues.procedure="x*"</literal>, then the algorithm is computing the solution
356                     of the linear program (phase 2).
357                 </para>
358             </listitem>
359         </itemizedlist>
360         <para>
361             It might happen that the output function requires additional arguments to be evaluated.
362             In this case, we can use the following feature.
363             The function <literal>outfun</literal> can also be the list <literal>(outf,a1,a2,...)</literal>.
364             In this case <literal>outf</literal>, the first element in the list, must have the header:
365         </para>
366         <programlisting role="example"><![CDATA[ 
367        stop = outf ( xopt , optimValues , state , a1 , a2 , ... )
368 ]]></programlisting>
369         <para>
370             where the input arguments <literal>a1, a2, ...</literal>
371             will be automatically be appended as parameters to the call.
372         </para>
373     </refsection>
374     <refsection>
375         <title>Stopping rule</title>
376         <para>
377             The stopping rule is based on the number of iterations,
378             the relative tolerance on the function value, the duality gap and 
379             the user's output function.
380         </para>
381         <para>
382             In both the phase 1 and phase 2 of the algorithm, we check the
383             duality gap and the boolean:
384         </para>
385         <programlisting role="example"><![CDATA[ 
386 dualgap > 1.e5 * dualgapmin
387 ]]></programlisting>
388         <para>
389             where <literal>dualgap</literal> is the absolute value of the duality gap, and 
390             <literal>dualgapmin</literal> is the minimum absolute duality gap measured during 
391             this step of the algorithm.
392             The duality gap is computed from 
393         </para>
394         <programlisting role="example"><![CDATA[ 
395 dualgap = abs(yopt'*beq - fopt)
396 ]]></programlisting>
397         <para>
398             where <literal>yopt</literal> is the Lagrange multiplier, <literal>beq</literal>
399             is the right hand side of the inequality constraints and <literal>fopt</literal>
400             is the minimum function value attained in the current phase.
401         </para>
402         <para>
403             During the second phase of the algorithm (i.e. once <literal>x0</literal> is
404             determined), the termination condition for the function value is based on the boolean:
405         </para>
406         <programlisting role="example"><![CDATA[ 
407     abs(fprev-fopt)<=rtolf*abs(fprev)
408 ]]></programlisting>
409         <para>
410             where <literal>fprev</literal> is the previous function value,
411             <literal>fopt</literal> is the current function value and
412             <literal>rtolf</literal> is the relative tolerance on the function value.
413         </para>
414     </refsection>
415     <refsection>
416         <title>Implementation notes</title>
417         <para>
418             The implementation is based on the primal affine scaling algorithm, as
419             discovered by Dikin in 1967, and then re-discovered by Barnes and Vanderbei et al in 1986.
420         </para>
421         <para>
422             If the scaling factor <literal>gam</literal> is closer to 1 (say <literal>gam=0.99</literal>,
423             for example), then the number of iterations may be lower.
424             Tsuchiya and Muramatsu proved that if an optimal solution exists, then, for any
425             <literal>gam</literal> lower than 2/3, the sequence converges to a point in the interior point of the optimal face.
426             Dikin proved convergence with <literal>gam=1/2</literal>.
427             Mascarenhas found two examples where the parameter <literal>gam=0.999</literal> lets the
428             algorithm converge to a vertex which is not optimal, if the initial guess is chosen
429             appropriately.
430         </para>
431     </refsection>
432     <refsection>
433         <title>Example #1</title>
434         <para>
435             In the following example, we solve a linear optimization problem with 2 linear
436             equality constraints and 3 unknowns.
437             The linear optimization problem is
438         </para>
439         <para>
440             <latex>
441                 \begin{eqnarray}
442                 \begin{array}{l}
443                 \textrm{minimize } -x_1 -x_2\\
444                 x_1 - x_2 = 0\\
445                 x_1 + x_2 + x_3 = 2\\
446                 x \geq 0
447                 \end{array}
448                 \end{eqnarray}
449             </latex>
450         </para>
451         <para>
452             The following script solves the problem.
453         </para>
454         <programlisting role="example"><![CDATA[ 
455     Aeq = [
456     1 -1 0
457     1  1 1
458     ];
459     beq = [0;2];
460     c = [-1;-1;0];
461     x0 = [0.1;0.1;1.8];
462     [xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c)
463     xstar=[1 1 0]'
464 ]]></programlisting>
465         <para>
466             The previous script produces the following output.
467         </para>
468         <screen><![CDATA[ 
469 -->[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c)
470  yopt  =
471    ineqlin: [0x0 constant]
472    eqlin: [2x1 constant]
473    lower: [3x1 constant]
474    upper: [3x1 constant]
475  iter  =
476     68.  
477  exitflag  =
478     1.  
479  fopt  =
480   - 1.9999898  
481  xopt  =
482     0.9999949  
483     0.9999949  
484     0.0000102  
485 ]]></screen>
486         <para>
487             We can explore the Lagrange multipliers by detailing
488             each field of the yopt structure.
489         </para>
490         <screen><![CDATA[ 
491 -->yopt.ineqlin
492  ans  =
493      []
494 -->yopt.eqlin
495  ans  =
496   - 6.483D-17  
497     1.         
498 -->yopt.lower
499  ans  =
500   - 2.070D-10  
501   - 2.070D-10  
502     1.         
503 -->yopt.upper
504  ans  =
505     0.  
506     0.  
507     0.  
508 ]]></screen>
509         <para>
510             We can as well give the initial guess x0, as in the following script.
511         </para>
512         <programlisting role="example"><![CDATA[ 
513     Aeq = [
514     1 -1 0
515     1  1 1
516     ];
517     beq = [0;2];
518     c = [-1;-1;0];
519     x0 = [0.1;0.1;1.8];
520     [xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,x0)
521 ]]></programlisting>
522         <para>
523             In the case where we need more precision, we can reduce the
524             relative tolerance on the function value.
525             In general, reducing the tolerance increases the number of iterations.
526         </para>
527         <screen><![CDATA[ 
528 -->[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,[],1.e-5);
529 -->disp([fopt iter])
530   - 1.9999898    68.  
531 -->[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,[],1.e-7);
532 -->disp([fopt iter])
533   - 1.9999998    74.  
534 -->[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,[],1.e-9);
535 -->disp([fopt iter])
536   - 2.    78.  
537 ]]></screen>
538     </refsection>
539     <refsection>
540         <title>Example #2</title>
541         <para>
542             In the following example, we solve a linear optimization problem with 10 random linear
543             equality constraints and 20 unknowns.
544             The initial guess is chosen at random in the [0,1]^p range.
545         </para>
546         <programlisting role="example"><![CDATA[ 
547 n=10;
548 p=20;
549 Aeq=rand(n,p);
550 c=rand(p,1);
551 x0=rand(p,1);
552 beq=Aeq*x0;
553 xopt=karmarkar(Aeq,beq,c,x0);
554 // Check constraints
555 norm(Aeq*xopt-beq)
556 ]]></programlisting>
557     </refsection>
558     <refsection>
559         <title>Inequality constraints</title>
560         <para>
561             Consider the following linear program with inequality constraints.
562         </para>
563         <para>
564             <latex>
565                 \begin{eqnarray}
566                 \begin{array}{l}
567                 \textrm{minimize } - 20 x_1 - 24 x_2\\
568                 3 x_1 + 6 x_2 \leq 60 \\
569                 4 x_1 + 2 x_2 \leq 32 \\
570                 \end{array}
571                 \end{eqnarray}
572             </latex>
573         </para>
574         <programlisting role="example"><![CDATA[ 
575 c = [-20 -24]';
576 A = [
577 3 6
578 4 2
579 ];
580 b = [60 32]';
581 xopt=karmarkar([],[],c,[],[],[],[],[],A,b)
582 ]]></programlisting>
583         <para>
584             The previous script produces the following output.
585         </para>
586         <screen><![CDATA[ 
587 -->xopt=karmarkar([],[],c,[],[],[],[],[],A,b)
588  xopt  =
589     3.9999125  
590     7.9999912  
591 ]]></screen>
592     </refsection>
593     <refsection>
594         <title>With bounds</title>
595         <para>
596             Consider the following linear optimization problem.
597             The problem is used in the Scilab port of the Lipsol 
598             toolbox by Rubio Scola (example0.sce).
599             The original Lipsol toolbox was created by Yin Zhang.
600         </para>
601         <para>
602             <latex>
603                 \begin{eqnarray}
604                 \begin{array}{l}
605                 \textrm{minimize } 2 x_1 + 5 x_2 - 2.5 x_3\\
606                 x_1 +   S4 x_3 \leq 5 \\
607                 E2 x1 - x2 - x3 \leq 0 \\
608                 -2 \leq x_1 \leq 2 \\
609                 1 \leq x_2 \\
610                 0 \leq x_3 \leq 3 \\
611                 \end{array}
612                 \end{eqnarray}
613             </latex>
614         </para>
615         <para>
616             where <literal>S4 = sin(pi/4)/4</literal> and <literal>E2 = exp(2)</literal>.
617         </para>
618         <programlisting role="example"><![CDATA[ 
619 c = [ 2; 5; -2.5];
620 S4 = sin(%pi/4)/4;
621 E2 = exp(2);
622 A = [
623 1  0 S4
624 E2 -1 -1
625 ];
626 b = [ 5; 0];
627 lb = [ -2; 1   ; 0 ];
628 ub = [  2; %inf; 3 ];
629 xstar = [-2;1;3];
630 [xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb,ub)
631 ]]></programlisting>
632         <para>
633             The previous script produces the following output.
634         </para>
635         <screen><![CDATA[ 
636 -->[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb,ub)
637     yopt  =
638     ineqlin: [2x1 constant]
639     eqlin: [0x0 constant]
640     lower: [3x1 constant]
641     upper: [3x1 constant]
642     iter  =
643     76.
644     exitflag  =
645     1.
646     fopt  =
647     - 6.4999482
648     xopt  =
649     - 1.9999914
650     1.0000035
651     2.9999931
652 ]]></screen>
653     </refsection>
654     <refsection>
655         <title>Configuring an output function</title>
656         <para>
657             It may be useful to configure a callback, so that we can customize the
658             printed messages or create a plot.
659             Consider the following linear optimization problem, which is
660             presented on Wikipedia 
661             in <ulink url="http://en.wikipedia.org/wiki/Karmarkar's_algorithm">Karmarkar's algorithm</ulink>.
662         </para>
663         <para>
664             <latex>
665                 \begin{eqnarray}
666                 \begin{array}{l}
667                 \textrm{minimize } - x_1 - x_2\\
668                 2 w x_1 + x_2 \leq 1+w^2, \quad w=0.0, 0.1, 0.2, ..., 1.0\\
669                 x_1,x_2 \geq 0
670                 \end{array}
671                 \end{eqnarray}
672             </latex>
673         </para>
674         <para>
675             The following output function plots the current point
676             and prints the iteration number, the value of the objective
677             function.
678         </para>
679         <programlisting role="example"><![CDATA[ 
680 function stop = myoutputfunction ( xopt , optimValues , state )
681     localmsg = gettext("Iter#%3.0f: %s (%s), f=%10.3e, x=[%s], gap=%10.3e\n")
682     xstr = strcat(msprintf("%10.3e\n",xopt)'," ")
683     mprintf(localmsg,optimValues.iteration,state,optimValues.procedure,..
684         optimValues.fval,xstr,optimValues.dualgap)
685     plot(xopt(1),xopt(2),"bo")
686     stop = %f
687 endfunction
688 ]]></programlisting>
689         <para>
690             The following script defines the optimization problem and
691             runs the optimization.
692         </para>
693         <programlisting role="example"><![CDATA[ 
694   n=11;
695   A = [2*linspace(0,1,n)',ones(n,1)];
696   b = 1 + linspace(0,1,n)'.^2;
697   c=[-1;-1];
698   // Plot the constraints
699   scf();
700   for i = 1 : n
701     plot(linspace(0,1,100),b(i)-A(i,1)*linspace(0,1,100),"b-")
702   end
703   // Run the optimization
704   xopt=karmarkar([],[],c,[],[],[],[],myoutputfunction,A,b);
705   // Plot the starting and ending points
706   plot(xopt(1),xopt(2),"k*")
707   ]]></programlisting>
708         <scilab:image>
709             function stop = myoutputfunction ( xopt , optimValues , state )
710             plot(xopt(1),xopt(2),"bo")
711             stop = %f
712             endfunction
713             
714             n=11;
715             A = [2*linspace(0,1,n)',ones(n,1)];
716             b = 1 + linspace(0,1,n)'.^2;
717             c=[-1;-1];
718             for i = 1 : n
719             plot(linspace(0,1,100),b(i)-A(i,1)*linspace(0,1,100),"b-")
720             end
721             xopt=karmarkar([],[],c,[],[],[],[],myoutputfunction,A,b);
722             plot(xopt(1),xopt(2),"k*")
723         </scilab:image>
724         
725         <para>
726             The previous script produces the following output and creates a graphics.
727         </para>
728         <screen><![CDATA[ 
729 -->xopt=karmarkar([],[],c,[],[],[],[],myoutputfunction,A,b);
730 Iter#  0: init (x0), f=1.000e+000, x=[0.000e+000 0.000e+000], gap=Inf
731 Iter#  0: init (x0), f=1.000e+000, x=[0.000e+000 0.000e+000], gap=Inf
732 Iter#  1: init (x0), f=5.000e-001, x=[2.201e-001 -4.313e-002], gap=3.676e-001
733 Iter#  2: init (x0), f=2.500e-001, x=[3.283e-001 -6.512e-002], gap=2.140e-001
734 Iter#  3: init (x0), f=1.250e-001, x=[3.822e-001 -7.634e-002], gap=1.161e-001
735 Iter#  4: init (x0), f=6.250e-002, x=[4.090e-001 -8.202e-002], gap=6.033e-002
736 Iter#  5: init (x0), f=3.125e-002, x=[4.225e-001 -8.488e-002], gap=3.072e-002
737 [...]
738 Iter# 50: init (x0), f=8.882e-016, x=[4.359e-001 -8.775e-002], gap=8.882e-016
739 Iter# 51: init (x0), f=4.441e-016, x=[4.359e-001 -8.775e-002], gap=4.441e-016
740 Iter# 52: init (x0), f=2.220e-016, x=[4.359e-001 -8.775e-002], gap=2.220e-016
741 Iter# 52: init (x*), f=-3.481e-001, x=[4.359e-001 -8.775e-002], gap=Inf
742 Iter# 52: iter (x*), f=-3.481e-001, x=[4.359e-001 -8.775e-002], gap=Inf
743 Iter# 53: iter (x*), f=-7.927e-001, x=[5.249e-001 2.678e-001], gap=5.098e-001
744 [...]
745 Iter# 65: iter (x*), f=-1.250e+000, x=[5.005e-001 7.494e-001], gap=1.258e-004
746 Iter# 66: iter (x*), f=-1.250e+000, x=[5.005e-001 7.494e-001], gap=5.941e-005
747 Iter# 67: iter (x*), f=-1.250e+000, x=[5.005e-001 7.495e-001], gap=2.882e-005
748 Iter# 68: iter (x*), f=-1.250e+000, x=[5.005e-001 7.495e-001], gap=1.418e-005
749 Iter# 69: done (x*), f=-1.250e+000, x=[5.005e-001 7.495e-001], gap=7.035e-006
750  xopt  =
751     0.5005127  
752     0.7494803  
753 ]]></screen>
754     </refsection>
755     <refsection>
756         <title>Configuring the scaling factor</title>
757         <para>
758             In some cases, the default value of the scaling factor <literal>gam</literal>
759             does not make the algorithm convergent. To make it converge, we may reduce the
760             value of <literal>gam</literal>. The next script shows two results: First with
761             <literal>gam=0.5</literal> and second with <literal>gam=0.3</literal>.
762         </para>
763         <programlisting role="example"><![CDATA[
764     A = [-0.1548,-0.0909,-0.0014,-0.0001;0.0989,-0.0884,0.0004,0];
765     B = [0.1966354;0.2167484];
766     C = [0.2056;0.0908;0.0012;0];
767     lb = [0;0;0;0];
768     //Test with default value (gam=0.5)
769     [xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],C,[],[],[],[],[],A,B,lb)
770     //Test reducing gam (gam=0.3)
771     gam=0.3;
772     [xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],C,[],[],gam,[],[],A,B,lb)
773 ]]></programlisting>
774         <para>
775             The previous script produces the following output.
776         </para>
777         <screen><![CDATA[
778 -->[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],C,[],[],[],[],[],A,B,lb)
779     yopt  =
780
781     ineqlin: [2x1 constant]
782     eqlin: [0x0 constant]
783     lower: [4x1 constant]
784     upper: [4x1 constant]
785     iter  =
786     94.
787     exitflag  =
788     - 2.
789     fopt  =
790     - 0.0003129
791     xopt  =
792     - 0.0010041
793     - 0.0011719
794     0.0000007
795     0.7218612
796
797 -->[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],C,[],[],gam,[],[],A,B,lb)
798     yopt  =
799
800     ineqlin: [2x1 constant]
801     eqlin: [0x0 constant]
802     lower: [4x1 constant]
803     upper: [4x1 constant]
804     iter  =
805     172.
806     exitflag  =
807     1.
808     fopt  =
809     0.0000005
810     xopt  =
811     0.0000016
812     0.0000015
813     2.811D-08
814     0.7349402
815 ]]></screen>
816     </refsection>
817     <refsection>
818         <title>Infeasible problem</title>
819         <para>
820             Consider the following linear optimization problem.
821             It is extracted from "Linear Programming in Matlab"
822             Ferris, Mangasarian, Wright, 2008, Chapter 3, "The Simplex Method", Exercise 3-4-2 1.
823         </para>
824         <programlisting role="example"><![CDATA[ 
825     // An infeasible problem.
826     // Minimize -3 x1 + x2
827     //  - x1 -   x2 >= -2
828     //  2 x1 + 2 x2 >= 10
829     // x >= 0
830     c = [-3;1];
831     A=[
832     -1 -1
833     2  2
834     ];
835     A=-A;
836     b=[-2;10];
837     b=-b;
838     lb=[0;0];
839     [xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb)
840 ]]></programlisting>
841         <para>
842             The previous script produces the following output.
843         </para>
844         <screen><![CDATA[ 
845 -->[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb)
846     yopt  =
847
848     ineqlin: [0x0 constant]
849     eqlin: [0x0 constant]
850     lower: [0x0 constant]
851     upper: [0x0 constant]
852     iter  =
853     40.
854     exitflag  =
855     - 1.
856     fopt  =
857     []
858     xopt  =
859     []
860 ]]></screen>
861     </refsection>
862     <refsection>
863         <title>Unbounded problem</title>
864         <para>
865             Consider the following linear optimization problem.
866             It is extracted from "Linear and Nonlinear Optimization"
867             Griva, Nash, Sofer, 2009, Chapter 5, "The Simplex Method", Example 5.3.
868         </para>
869         <programlisting role="example"><![CDATA[ 
870   // An unbounded problem.
871   // Minimize -x1 - 2 x2
872   //  -1 x1 +   x2 <= 2
873   //  -2 x1 +   x2 <= 1
874   // x >= 0
875   c = [-1;-2];
876   A=[
877   -1  1
878   -2  1
879   ];
880   b=[2;1];
881   lb=[0;0];
882   [xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb)
883 ]]></programlisting>
884         <para>
885             The previous script produces the following output. 
886             Notice that the function produces <literal>exitflag=-2</literal>, 
887             which indicates that the algorithm detects that the duality gap has increased much
888             more than expected. 
889             This is the sign for a failure of the algorithm to find an optimal point.
890         </para>
891         <screen><![CDATA[ 
892 -->[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb)
893     yopt  =
894     ineqlin: [2x1 constant]
895     eqlin: [0x0 constant]
896     lower: [2x1 constant]
897     upper: [2x1 constant]
898     iter  =
899     59.
900     exitflag  =
901     - 2.
902     fopt  =
903     - 45374652.
904     xopt  =
905     15124883.
906     15124885.
907 ]]></screen>
908     </refsection>
909     <refsection>
910         <title>References</title>
911         <para>
912             "Iterative solution of problems of linear and quadratic programming",
913             Dikin,
914             Doklady Akademii Nausk SSSR, Vol. 174, pp. 747-748, 1967
915         </para>
916         <para>
917             "A New Polynomial Time Algorithm for Linear Programming",
918             Narendra Karmarkar, Combinatorica, Vol 4, nr. 4, p. 373–395, 1984.
919         </para>
920         <para>
921             "A variation on Karmarkar’s algorithm for solving linear programming problems,
922             Earl R. Barnes, Mathematical Programming, Volume 36, Number 2, 174-182, 1986.
923         </para>
924         <para>
925             "A modification of karmarkar's linear programming algorithm",
926             Robert J. Vanderbei, Marc S. Meketon and Barry A. Freedman,
927             Algorithmica, Volume 1, Numbers 1-4, 395-407, 1986.
928         </para>
929         <para>
930             "Practical Optimization: Algorithms and Engineering Applications",
931             Andreas Antoniou, Wu-Sheng Lu, Springer, 2007,
932             Chapter 12, "Linear Programming Part II: Interior Point Methods".
933         </para>
934         <para>
935             "Global Convergence of a Long-Step Affine Scaling Algorithm for Degenerate Linear Programming Problems",
936             Takashi Tsuchiya and Masakazu Muramatsu,
937             SIAM J. Optim. Volume 5, Issue 3, pp. 525-551 (August 1995)
938         </para>
939         <para>
940             "The convergence of dual variables",
941             Dikin,
942             Tech. report, Siberian Energy Institute, Russia, 1991
943         </para>
944         <para>
945             "The Affine Scaling Algorithm Fails for Stepsize 0.999",
946             Walter F. Mascarenhas, SIAM J. Optim. Volume 7, Issue 1, pp. 34-46 (1997)
947         </para>
948         <para>
949             "A Primal-Dual Exterior Point Algorithm For Linear Programming Problems"
950             Nikolaos Samaras, Angelo Sifaleras, Charalampos Triantafyllidis
951             Yugoslav Journal of Operations Research
952             Vol 19 (2009), Number 1, 123-132
953         </para>
954     </refsection>
955 </refentry>