View difference between Paste ID: QEaiTpdx and bmpPJ83U
SHOW: | | - or go back to the newest paste.
1
function [x_opt,y_opt,W]=PSO(epsilon,Nmax,xmin,xmax)
2
3
W=0;
4
N=2; %wymiar problemu, tu dla dwóch zmiennych
5
C=20;
6
w=0.8; % wartosc powinna malec w czasie dzialania algorytmu
7
c1=1;
8
c2=1;
9
10
R.x=zeros(N,C);
11
R.v=zeros(N,C);
12
R.y=zeros(1,C);
13
14
for i=1:C
15
    R.x(:,i)=(xmax-xmin).*rand(N,1)+xmin;
16
    R.v(:,i)=0.02*(xmax-xmin).*rand(N,1)-0.01*(xmax-xmin);
17
    [R.y(i),W] = fit_fun(R.x(:,i),W);
18
end
19
20
x_opt_i=R.x;
21
y_opt_i=R.y;
22
[y_opt,numer]=min(R.y);
23
x_opt=R.x(:,numer);   %najlepsze rozwiązanie to cząstka o numerze ... 
24
25
while true
26
    for i=1:C
27
        R.v(:,i)=w*R.v(:,i)+c1*rand*(x_opt-R.x(:,i))+...
28
            c2*rand*(x_opt_i(:,i)-R.x(:,i));
29
        R.x(:,i)=R.x(:,i)+R.v(:,i);
30
        [R.y(i),W]=fit_fun(R.x(:,i),W);
31
        if R.y(i)<y_opt
32
            y_opt=R.y(i);
33
            x_opt=R.x(:,i);
34
            if y_opt<epsilon
35
                return;
36
            end
37
        end
38
        if R.y(i) <y_opt_i(i)
39
            y_opt_i(i)=R.y(i);
40
            x_opt_i(:,i)=R.x(:,i);
41
        end
42
        if W>Nmax
43
            return;
44
        end
45
    end
46
end
47
48
49
function [y,w]=fit_fun(x,w)
50
51
w=w+1;
52
%y=x(1)^2+x(2)^2;
53
m1=400;
54
m2=50;
55
k1=x(1);
56
k2=2e5;
57
b=x(2);
58
g=9.81;
59
x2_0=-(m1+m2)*g/k2;
60
x1_0=-m1*g/k1+x2_0;
61
T=5;
62
dt=0.01;
63
A=0.06;
64
f=2*pi*10;
65
66
assignin('base','m1',m1);
67
assignin('base','m2',m2);
68
assignin('base','k1',k1);
69
assignin('base','k2',k2);
70
assignin('base','b',b);
71
assignin('base','g',g);
72
assignin('base','x1_0',x1_0);
73
assignin('base','x2_0',x2_0);
74
assignin('base','T',T);
75
assignin('base','dt',dt);
76
assignin('base','a',A);
77
assignin('base','f',f);
78
sim('model.mdl');
79
plot(t,x1);
80-
y=0; % obliczenie wartosci funkcji celu na podstawie x1 ZROBIĆ !!!!! 
80+
figure;
81
82
[l,m]=butter(12,0.16,'low');
83
x1_f=filter(l,m,x1);
84
plot(t,x1_f);
85
86-
%wykres sygnału filtrowanego
86+
y=max(abs(x1_f));
87
y=y+1e6*p;
88
p=max(0,-k1)^2+max(0,k1-2e4)^2+ ...
89
    max(0,b)^2+max(0,b-3e3)^2+ ...
90-
figure
90+
91
plot(t,x1);
92-
y=mean(max(x1_f)+min(x1_f));
92+
93
94
95
96
error('To by było na tyle');