-
Notifications
You must be signed in to change notification settings - Fork 5
/
pwp_approximate.m
86 lines (73 loc) · 1.87 KB
/
pwp_approximate.m
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
function PWP = pwp_approximate(NLFcn)
% Usage: PWP = pwp_approximate(NLFcn)
% Example:
% NLFcn.Handle = @nlfcn;
% NLFcn.Resolution = Resolution;
% NLFcn.Domain = Domain;
% NLFcn.Order = nl;
%
% PWP = pwp_approximate(NLFcn);
% coeff = PWP.coeff;
% X = PWP.X;
% NLFun = PWP.NLFun;
% PWPFun = PWP.PWPFun;
%
% figure(100);
% line(X,NLFun,'Color',[0 0 0])
% line(X,PWPFun,'Color',0.7*[1 1 1],'LineWidth',2,'LineStyle','--');
% xlabel('$x$');
% ylabel('$f(x)$');
% legend('$f(x)$','$\hat f(x)$',2);
% set(gca,'Box','on','XTick',Domain,'XLim',[Xmin Xmax],'XGrid','on');
Domain = NLFcn.Domain;
Resolution = NLFcn.Resolution;
nl = NLFcn.Order;
xstar = NLFcn.xstar;
NR = length(Domain)-1;
X = sort(union(linspace(Domain(1),Domain(end),Resolution),Domain));
F = [];
x = sdpvar;
[p,c,v] = polynomial(x,nl);
P = sdpvar(1,length(v));
Funx = char(sdisplay(P*v));
FunX = strrep(Funx,'x','X(In)');
FunX = strrep(FunX,'^','.^');
% Approximation Error
for i = 1:NR,
In = find(X>=Domain(i)&X<=Domain(i+1));
NLFun(In) = feval(NLFcn.Handle,X(In));
coeff{i} = sdpvar(1,length(v));
P = coeff{i};
PWPFun(In) = eval(FunX);
end
Err = PWPFun-NLFun;
Obj = sdpvar(1);
F = [F, [Obj Err;Err' eye(length(Err))]>0];
% Continuity
for i=1:NR-1,
x = Domain(i+1);
P = coeff{i};
PWPFuni = eval(Funx);
P = coeff{i+1};
PWPFunip1 = eval(Funx);
F = [F, PWPFuni==PWPFunip1];
end
Istar = min(find(Domain>=xstar))-1;
x = xstar;
P = coeff{Istar};
PWPFuni = eval(Funx);
NLFunc = feval(NLFcn.Handle,x);
F = [F, PWPFuni==NLFunc];
sol = solvesdp(F,Obj);
for i = 1:NR,
In = find(X>=Domain(i)&X<=Domain(i+1));
NLFun(In) = feval(NLFcn.Handle,X(In));
coeff{i} = double(coeff{i});
P = coeff{i};
PWPFun(In) = eval(FunX);
end
PWP.coeff = coeff;
PWP.X = X;
PWP.NLFun = NLFun;
PWP.PWPFun = PWPFun;
PWP.NR = NR;