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
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
//
function [fpen,gpen,ind]=pencost(x,ind,fncts,ne,nc,cpen);
// pencost=template of external function for penalized problem:
//
// min f(x)
// g_i(x) = 0 i=1:ne
// g_i(x) <= 0 i=ne+1:nc (0<=ne<=nc)
//
//function pencost defines a penalized cost function and its gradient
//which is passed to scilab's optim function.
//
//A penalized problem is defined through the scilab function fncts
//which defines the constraints g_1,..,g_nc and the cost function f.
//
//fncts = scilab external with calling sequence:
//[gisandf,theirgrads,indic]=fncts(x,indic)
//which returns
// in row vector gisandf [g_1(x),...,g_nc(x),f(x)]
//
// in column 1 of matrix theirgrads grad(g_1)_x
//..........
// in column i of matrix theirgrads grad(g_i)_x
//...........
// in column nc of matrix theirgrads grad(g_nc)_x
// in column nc+1 of matrix theirgrads grad(g_nc)_x
//
// ne = number of equality constraints
// nc = total number of constraints
//
//EXAMPLE:
//
//minimize f(x)=0.5*(x(1)^2+x(2)^2)
// g_1(x)= x(1)+x(2)-1 <=0
// g_2(x)=- x(1) <= 0
//
// no equality constraint (ne=0)
// two inequality constraints (-> nc=2 constraints)
//
// 1-define your problem through a function (usually in a .sci file):
//deff('[gsf,grds,indic]=myproblem(x,indic)',...
// 'gsf=[x(1)+x(2)-1,-x(1),norm(x)^2];...
// grds=[1,-1,x(1);...
// 1,0,x(2)]')
//
// 2-call scilab's optim function (after loading pencost.sci (this file));
// exec('pencost.sci','c')
// ne=0;nc=2;x0=[4;7];
// penalizing constant cpen=100;
// [f,x]=optim(list(pencost,myproblem,ne,nc,cpen),x0)
//
// or (passing myproblem,ne,nc,cpen by context)
// fncts=myproblem;[f,x]=optim(pencost,x0)
//
[f,gradf,indic]=fncts(x,ind);
if indic < 0 then ind=indic, return, end;
if nc >ne then
for i=ne+1:nc, f(i)=max([0 f(i)]); end;
end;
fpen=f(nc+1) + 0.5*cpen*norm(f(1:nc))^2;
if ind==2 then return, end;
gpen=gradf(:,nc+1);
if ne >0 then
for i=1:ne, gpen=gpen + cpen*f(i)*gradf(:,i),end;
end;
if nc > ne
for i=ne+1:nc,
if f(i)>0 then gpen=gpen + cpen*f(i)*gradf(:,i),end;
end;
end;
endfunction
|
