1/12/2007

好的matlab的網站

http://matlabdb.mathematik.uni-stuttgart.de/index.jsp

吳子青: Radial Basis Function ANN 輻狀基底類神經網路

Radial Basis Function ANN 輻狀基底類神經網路

前言:
利用RBF ANN 作簡短的介紹類神經網路的原理,並以老師上課的例子以類神經網路作簡單的應用,在這裡利用ANN 作 “ 函數模擬 ” 。

類神經網路介紹:
類神經網路(Artificial Neural Networks, ANNs)或譯為人工神經網路,其主要的基本概念是嘗試著模仿人類的神經系統。

它是由很多非線性的運算單元(即:神經元 neuron)和位於這些運算單元間的眾多連結(links)所組成,而這些運算單元通常是以平行且分散的方式來進行運算以電腦的軟硬體來模擬生物神經網路的資訊處理系統。

整個ANN的聚集形式就如同人類的大腦一般,可透過樣本或資料的訓練來展現出學習(learn)、回想(recall)、歸納推演(generalize)的能力。

類神經網路在處理複雜的工作時
(1)不需要針對問題定義複雜的數學模式,
(2)不用去解任何微分方程、積分方程或其他的數學方程式,
(3)藉由學習來面對複雜的問題與不確定性的環境。































輻狀基底函數類神經網路介紹(RBF ANN)
或稱為半徑式類神經網路,特質主要在於摩擬大腦皮質軸突的局部調整功能
為基本前饋式類神經網路
其可視為在解決高維度空間的曲線調適問題

RBFANN 架構


RBFANN演算流程
1. 處理資料 (輸入項與目標輸出項)
2. 決定類神經網路架構
3. 決定初始中心點位置 (類神經元)
4. 利用RBFANN學習演算法修正(1)中心點位置 (2)連結權重 (3)福狀基底半徑
5. 驗證與測試階段


程式流程
Step 1 :
想要利用RBF ANN模擬 matlab裡面z=peaks(X,Y)的函數
z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ...
- 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ...
- 1/3*exp(-(x+1).^2 - y.^2)
可以看出 Peaks(X,Y) 這是一個非線性的方程式
正好用來測試類神經網路處理非線性的能力
取X,Y都介於-2和2間距0.2 共產生400筆資料各包含(X,Y, f(X,Y ))

Step 2
架構RBFANN 這裡選取輸入維度2, 中心點個數 100
以高斯函數為輻狀基底函數 演算次數 200次 學習速率0.01

Step3
將400筆資料依序丟入網路中進行訓練,並修正中心點位置 連結權重 福狀基底半徑 等資訊

Step4
製作演算流程動畫, 以3D圖上的黑點表示 input data , 以mesh 網格表示ANN模擬函數的結果, 若資料點和mesh網格重合則表示模擬效果優良

Step 5
達到停止條件或最大演算步數則停止, 結果為rmse 越小表示越好

Step 6
製作GUI介面
程式內容
clear all;
iter=100;
nc=100;
goal=0.01;
lr_w=0.01;
lr_c=0.01;
lr_s=0.01;

% p 為輸入資料點,N×K矩陣,N是輸入資料維度,K是資料點數
% t 為目標輸出值,1×K矩陣
% newcenter已選定的中心點,N×nc矩陣
% iter為指定演算代數
% goal為指定之網路誤差
% lr_w為調整權重之學習速率
% lr_c為調整中心點之學習速率
% lr_s為調整標準偏差? (sigma)之學習速率
% W為輸出層權重,nc×1矩陣
% yh為網路輸出值,1×K矩陣
% rmse為目標輸出值與網路輸出值之RMSE

%輸入函數
[X,Y] = meshgrid(-2:0.2:2);
Z = 3*peaks(X,Y);

[mm,nn]=size(Z);
PA=zeros(3,mm*nn);
for i=1:mm
for j=1:nn
PA(1,(i-1)*nn+j)=X(i,j);
PA(2,(i-1)*nn+j)=Y(i,j);
PA(3,(i-1)*nn+j)=Z(i,j);
end
end
clear i j ;

p=PA(1:2,:);
t=PA(3,:);

[nd,np]=size(p);
k=randperm(np);
newcenter=p(:,k(1:nc));

clear nd np nc k ;

nc=size(newcenter,2);
[nd,np]=size(p);
phii=zeros(1,nc);
Dic=[];
phi=[];
itimetemp=1;
for i=1:nc
for j=1:nc
Dic(j,i)=norm(newcenter(:,i)-newcenter(:,j));%計算center的距離
end
end
sigma=(max(max(Dic))/sqrt(nc));%計算基底函數的標準偏差σ(sigma)
for i=1:nc
for j=1:np
phi(j,i)=exp(-(norm(newcenter(:,i)-p(:,j))/sigma)^2);%計算center與input各點的距離
end
end

W=pinv(phi)*t';%初始權重W
newsigma=sigma*ones(nc,1);
for itime=1:iter
for j=1:np
for i=1:nc
Dip(i)=norm(p(:,j)-newcenter(:,i));%計算center與input各點的距離
phii(i)=exp((-1/(newsigma(i)^2))*(Dip(i)^2));
end
if(itime==iter)
finalphii(j,:)=phii;
end
yh(j)=phii*W;
e(j)=yh(j)-t(j);
for k=1:nc
newcenter(:,k)=newcenter(:,k)-lr_c*(e(j)*W(k)/(newsigma(k)^2))*phii(k)*(p(:,j)-newcenter(:,k));%更新基底函數的中心值c
newsigma(k)=newsigma(k)-lr_s*(e(j)*W(k)/(newsigma(k)^3))*...
phii(k)*(norm(p(:,j)-newcenter(:,k)))^2;%更新基底函數的標準偏差σ
end
W=W-lr_w*e(j)*phii';%更新權重向量W
end
err1=sqrt(mse(yh-t));
err2(itime)=mse(yh-t);
if(itime/100-itimetemp==0)
itimetemp=itimetemp+1;
itime
end
if(err1<=goal)
break;
end

for i=1:np
for m=1:nc
phii(m)=exp(-(norm(p(:,i)-newcenter(:,m))/newsigma(m))^2);
end
yh(i)=phii*W;
end
for i=1:mm
for j=1:nn
ZZ(i,j)=yh((i-1)*nn+j);
end
end
rmse=sqrt(mse(yh-t));
clear yh;
%figure1
figure(1)
plot3(PA(1,:),PA(2,:),PA(3,:),'.','markersize',5);
grid on
hold on
mesh(X,Y,ZZ);
hold off
M(itime)=getframe;
RMSE(itime)=rmse;
end %time end
movie(M,1,2);

figure(2)
plot(RMSE);
title('rmse value');
xlabel('epoch step');
ylabel('rmse')
clear Dic Dip W ans e err1 err2 finalphii goal itime itimetemp;
clear j k lr_c lr_s lr_w m nc nd newsigma np p phi phii sigma itime m n mm nn;
clear PA X Y Z ZZ newcenter;
clear itre i M iter t RMSE;

林書如: 人體受外力牽拉後引發之姿勢反應_肌肉活化情形_finale

請按標題連結至程式網頁.....

林翊展--養豬廢水處理設備的自動化設計--完整版

前言
前一陣子曾經協助老師進行一項雲林縣政府的計畫,該計畫內容為替雲林縣政府設計一座養豬廢水處理的示範場.在設計過程當中,有需要計算的地方,雖然用到的計算方法不難,但是希望利用matlab程式,只需要輸入所需的參數,即可以自動設計出符合標準的處理設備規格.

背景說明

養豬廢水處理

台灣的畜牧業以養豬為大宗,而且多是採企業化經營,規模大又集中,導致畜牧場所產生的廢水量相當驚人.目前台灣的總養豬頭數大約700萬頭,大部分都集中在台灣中南部,而一頭豬的糞尿廢水汙染力相當於人的六倍,可以折算成約4200萬人口.在過去,養豬廢水都是以直接排放的方式處理,對於台灣河川的破壞相當大,產生了許多”黑水河”.現在依環保署的規定,養豬規模超過20頭要強制建設廢水處理設備,而處理過的廢水也有品質的要求,如此一來,大大減少台灣河川所遭受的污染,而且在處理豬糞尿的過程中,還會產生許多具有附加價值的產物,如沼氣、有機堆肥等,達到保護環境與廢棄物再利用的雙重功效.

目前國內豬糞尿水的處理一般都採用三段式處理,過程主要包括固液分離、厭氣醱酵及好氣處理等步驟,其中再加入調整池、計量槽、初沉池、污泥回流及污泥處理等構成完整的養豬廢水處理系統,其流程如下圖:




而實際處理圖片則如下所示

1.固液分離




2.厭氣發酵




3.好氣發酵




經過三段式處理後的廢水即放流到河川之中

放流廢水




目的

程式設計的目標為輸入養豬頭數後,自動計算出廢水處理槽的規格參數,節省人力計算所需的時間.


名詞解釋
1.BOD:biochemical oxygen demand,生化需氧量.水中的雜質經由好氧菌分解所需要的氧氣量.由此數值可以看出廢水的生物可分解性,亦是廢水污染力的重要指標之一.BOD5則是指讓廢水暴露在氧氣充足的情況下,五天所耗費的氧氣量.

2.COD:chemical oxygen demand,化學需氧量.利用化學方法(強氧化劑),讓水中雜質氧化的耗氧量.廢水污染指標之一,且COD值必大於BOD值,而COD值通常被當作終極BOD值來作計算.

3.HRT:水力停留時間
HRT= reactor volume / effluent volume flux

4.BOD loading:反應槽對於進流廢水的負載,可以用BOD或COD表示,通常為經驗值,過多或過少都會對反應槽的效率產生影響.


設計原則
1.假設一頭肉豬100 kg,產生30L/day的廢水量
2.根據ASAE在1994公佈的STANDARDS 裡,每1000kg豬產生的BOD5約 3.1 kg/day.
3.反應槽的BOD loading假設為0.5 kg/m^3/day


程式內容


k=menu('請選擇計算方式','依廢水體積計算','依BOD負載計算');
pig=input('請輸入豬隻頭數 : ');
HRT=input('請輸入水力停留時間(hour) : ');
if k==1
q=pig*30*0.001;
V=q*(HRT/24);
disp(['反應槽體積= ',num2str(V),' 立方公尺']);
elseif k==2
V=pig*0.31/0.5*(HRT/24);
disp(['反應槽體積= ',num2str(V),' 立方公尺']);
else
break
end



執行結果
程式計算方式有兩種,第一是根據廢水量來計算,第二則是根據預估的反應槽loadind來做計算.兩種結果算出來後,取其較大值作為依據.

選擇圖示



根據廢水量計算



根據反應槽的loading來計算




根據以上的結果,養豬1000頭,水力停留時間6 hours的廢水處理槽,體積需要155立方公尺.


結果討論
1.兩種計算方式的數據會有不同,因為第一個根據廢水量的計算,只單就以物理性而言,槽體的大小是為了不讓水滿出來,但是第二項則是實驗的經驗值,參雜了細菌生長的好壞等生物性因素,所以兩者的參數不相同,必須取其大者.

2.水力停留時間6~8小時,是傳統反應槽的設計,如果有需要,例如loading改變,細菌生長狀況差異等,則會另行調整,以便達到處理效率.

3. 現行環保署畜牧業放流水標準
依據水污染防制法第七條第二項之規定,國內畜牧業適用之放流水標準,其水質項目及限制如下:
(1)非草食性動物(豬):生化需氧量 < 80 mg/L;化學需氧量 < 600 mg/L;懸浮固體 < 150 mg/L。
(2)草食性動物(牛):生化需氧量 < 80 mg/L;化學需氧量 < 450 mg/L;懸浮固體 < 150 mg/L。
(行政院環保署,2003)

4.以現實情況來說,廢水處理達不到標準通常受限於成本與土地大小,造成反應槽的大小不足,再加上生物性因素會比較難以掌握.

5.固液分離後的豬糞可以進行發酵做為有機堆肥,而厭氧反應過程中會產生沼氣提供發電,不過要先想辦法去除沼氣中的硫化氫,因為它會產生硫酸腐蝕發電機組.

6.此篇的程式與其他同學比較起來相對的簡單很多,因為生物反應槽大部分都是實驗參數,設計上比較難的地方在於參數的取用,用到計算上都只是簡單的計算,而現實上則隨時根據生物狀況來做調整,不過此篇可以提供同學們一點關於廢水處理的相關資訊.

1/11/2007

r95631022林冠宏(計算齒輪平均半徑,有效圈數,彈簧體積)

螺旋彈簧
緊密捲繞(closely coiled)的螺旋彈簧(helical springs)應力與變形方程式,直接導自如圖4-2(2)中所示圓桿扭轉的對應力方程式。圓桿長L,直徑為d,在兩端裝配長R的托架,並且在負荷P的作用下平衡。假設直桿彎曲成如圖4-2(b)中半徑為R的N圈螺旋線(helix)。繞成圈的桿在兩相等兒相反方向力的作用下平衡。
直桿中的剪應力源自值等於PR的扭矩。螺旋線中的主要應力也是扭轉剪應力。











由扭轉剪應力=Tr/J(J為極慣性矩),得圖4-2中的應力為


(a)

將螺旋直徑2R與線徑d間的比值為彈簧指數(spring index) c。即



(1)

在繞成螺旋線後,由橫向剪力的作用,剖面將有額外應力。經精確分析顯示,此應力在彈簧長的中央部份之值為1.23P/A。則



(b)



由於靜負荷P的作用,在彈簧長中央部份的線圈內側,其總剪應力τ可由(a)式與(b)是相加求得:

(2)


式中

(3)


將(1)式代入(2)式得另外兩個型式的應力方程式

(4)(5)


當(3)式代入(4)式而解得c,可得下列有用的方程式

(6)

將(1)式的c代入(6)式,得

(6a)


彈簧的撓度可由考慮由扭矩PR導致剖面間的互相旋轉而得。暫時假設圖4-2(b)中的元素ABCD據撓性,彈簧的其他部分為剛體。則φ=Tl/JG式,剖面CD對相鄰剖面AB的旋轉dφ 等於



此一微分長度彈簧得旋轉,導致位於相隔 的E點的總移動量為




當整個彈簧具彈性時,由扭矩產生的總撓度,可以由全長積分而得,則

(7)(8)(9)


在這些方程式中,彈簧的作用長度(active length)取為2πRNc,而鋼線極貫性矩為πd^4/32。方程式(8)及(9)由(1)式代入而得。
彈簧用鋼的G,平均值為11500000 psi ,若以SI單位表示。為79300MPa。
若Q為彈簧的無效端圈數,則彈簧總長為2πr(Nc+Q)。鋼線的剖面積為πd^2/4彈簧材料的體積V為兩項的乘積。

(10)

彈簧率(spring rate)k,或產生單位撓度所須之力的方程式可由將(7),(8),(9)式中的P以k,而δ以1取代而得


(11)


其他彈簧率k公式,可由考慮圖4-2(c)而得

(12)

題目:
由(幾)號剛線製成的螺旋壓縮彈簧,承受(多少lb)靜負荷形成1in撓度。及彈簧所設計承受最大剪應力將為(多少psi)。而無效端圈數Q等於(多少)。可以知道所需平均半徑R,有效圈數 及彈簧材料的體積之值。
分析:
由圖表1.1知











表4-1

由表4-1得線徑d
由(6a)式知R
由(1)式知C
由(9)式知δ
得Nc
由(10)式5 之V

% gear.m
% Design of a gear
ld=[3/8 0.375
5/16 0.3125
1/4 0.25
4 0.2253
5 0.207
6 0.192
7 0.177
8 0.162
9 0.1483
10 0.135
11 0.1205
12 0.1055
13 0.0915
14 0.08
15 0.072
16 0.0625
17 0.054
18 0.0475
19 0.041
20 0.0348
21 0.0317
22 0.0286
23 0.0258
24 0.023
25 0.0204];
while 1

No=input('幾號鋼線製成的螺旋壓縮彈簧(-)[預設值為8]:','s');
if isempty(No)
No='8';
end
number=str2num(No);
[N,I]=find(ld(:,1)==number);


if isempty(N) %找表中較近的值
if N>25
while isempty(N)
number=number-1;
[N,I]=find(ld(:,1)==number);
end
disp([' 表中較近的值為 ' num2str(number)]);
yn= upper(input('是否使用(Y/N) [Y] ','s'));
if yn ~= 'Y'
number=0;
[N,I]=find(ld(:,1)==number);
end
end


end


if isempty(N)
disp('並無此號數彈簧!');
else
break;
end
end

while 1
F=input('請輸入承受靜負荷形成1in撓度(lb)[預設值為100]:','s');
if isempty(F)
F='100';
end
P=str2num(F);


T=input('最大剪力(psi)[預設值為80000]:','s');
if isempty(T)
T='80000';
end
t=str2num(T);


Q=input('請輸入無效半徑[預設值為2]: ','s');
if isempty(Q)
Q='2';
end
q=str2num(Q);

d=ld(N,2);
R=d/2*((3.14*d^2*t/8/P)-0.615);
c=2*R/d;
G=11500000;
nc=1*R*G/4/P/c^4;
v=1/2*pi^2*d^2*R*(nc+q);

disp('齒輪資料:');
disp(['第' num2str(N) '號彈簧']);
disp(['承受靜負荷形成1in撓度:' num2str(F) 'lb']);
disp(['最大剪力(psi):' num2str(T) 'psi']);
disp(['無效半徑:' num2str(Q)]);
disp(['所需平均半徑R:' num2str(R) 'in']);
disp([' 有效端圈數N:' num2str(nc)]);
disp(['彈簧材料的體積:' num2str(v) 'in3.']);

A=upper(input('繼續執行嗎?(Y/N) [Y] ','s'));
if A~='Y'
break;
end
end
結果:


請輸入承受靜負荷形成1in撓度(lb)[預設值為100]:300
最大剪力(psi)[預設值為80000]:30000
請輸入無效半徑[預設值為2]: 2
齒輪資料:
第5號彈簧
承受靜負荷形成1in撓度:300lb
最大剪力(psi):30000psi
無效半徑:2
所需平均半徑R:0.11042in
有效端圈數N:816.92
彈簧材料的體積:19.1199in3.
繼續執行嗎?(Y/N) [Y]

1/04/2007

熱傳模組建立(續)非穩態的建立

這次原來的帳號似乎官方好像修好了,於是又改回在原gmail的申請部落格寫
(請點標題連結)