Appendix A. MATLAB Scripts

Script 1: graphing functions u1 and u2 in design space in Example 2.9.
clear all;
[q1,q2] = meshgrid(0:1:100,0:1:100);
u1 = -q1.ˆ 2+(198-q2).∗q1;
[C,h] = contour(q1,q2,u1,3);
set(h,'ShowText','on','TextStep',get(h,'LevelStep')∗2)
colormap cool
hold on;
u2 = -q2.ˆ 12+(198-q1).∗q2;
[C,h] = contour(q1,q2,u2,3);
set(h,'ShowText','on','TextStep',get(h,'LevelStep')∗2)
colormap cool.
hold off;
Script 2: graphing functions u1 and u2 in criterion space in Example 2.9
%Cournot example; plot objective domain
clear all; format long
NN = 10,000;  %No. of random points for plotting
%----------------------------------------------
A = 1; %counting feasible points
B = 1; %counting infeasible points
Smax = 0;
q1o = 0;
q2o = 0;
for ii = 1:NN
    q1 = rand∗100;
    q2 = rand∗100;
    u1 = -q1 ˆ 2+(198-q2)∗q1;
    u2 = -q2 ˆ 2+(198-q1)∗q2;
    if q1>0 & q2>0;
        feasi(:,A) = [u1,u2];
        if Smax < feasi(1,A);
            Smax = feasi(1,A);
            q1o = q1; q2o = q2;
        end
        A = A+1;
    else
        infeasi(:,B) = [u1,u2];
        B=B+1;
    end
end
figure (1)
plot(feasi(1,:),feasi(2,:),'bs','MarkerSize',3)
plot(infeasi(1,:),infeasi(2,:),'rs','MarkerSize',3)
xlabel('u1′); ylabel('u2′)
axis equal
Script 3: solving the unconstrained problem of the cantilever beam example in Section 2.6.1.1
%Unit: mm,N,MPa
clear all
%------------- Parameters that can be adjusted
k_w = 0.5; k_z = 0.5;        %%%%% Weighting between attributes
r_w = 0.0001;  r_z = −0.0001;  %%%%% Model risk; r > 0 risk adverse; r < 0 risk adverse; r∼0 risk neutral
b_u = 60; b_l = 30;         %%%%% Upper and lower bounds of base
h_u = 80; h_l = 40;         %%%%% Upper and lower bounds of height
interval = 0.1;           %%%%% Solution search interval size
%------------- Define properties
l=200;         P=1000;    E=69,000;    YS=27.6;
rou = 2700∗10 ˆ (-9);  %Kg/mm3
%------------------------- An example; special case
b=30;h=40;I=b∗hˆ3/12; z=∗lˆ3/(3∗E∗I);M=P∗l; stress=M∗(h/2)/I;w=rou∗b∗h∗l∗9.8;
%-----------------------------------
wl=rou∗b_u∗h_u∗l∗9.8;          w_u = rou∗b_l∗h_l∗l∗9.8;
z_l=P∗lˆ3/(3∗E∗(b_l∗h_lˆ3/12));      z_u= P∗lˆ3/(3∗E∗(b_u∗h_uˆ3/12));
 
x_w = w_l: 0.01:w_u;           x_z = z_l: 0.001:z_u;
s_w=(x_w-w_l)./(w_u-w_l);        s_z=(x_z-z_l)./(z_u-z_l);
u_w=(1-exp(-r_w.∗s_w))./(1-exp(-r_w));  u_z=(1-exp(-r_z.∗s_z))./(1-exp(-r_z));
 
figure (1);  % plot utility curves
subplot(1,2,1); plot(x_w,u_w); xlabel('Self weight'), ylabel('Utility')
subplot(1,2,2); plot(x_z,u_z); xlabel('Displacement'), ylabel('Utility')
%------------------------------------------ Build MAU
B = b_l:interval:b_u;
H = h_l:interval:h_u;
W = rou∗B.'∗H∗l∗9.8;
Z = P.∗l.ˆ3./(3∗E∗(B.'∗H.ˆ3/12));
 
S_w=(W-w_l)./(w_u-w_l);
S_z=(Z-z_l)./(z_u-z_l);
 
U_w=(1-exp(-r_w.∗S_w))./(1-exp(-r_w));
U_z=(1-exp(-r_w.∗S_z))./(1-exp(-r_w));
K=(1-k_w-k_z)/(k_w∗k_z);
U_MAU = k_w∗U_w + k_z∗U_z + K∗k_w∗k_z∗U_w.∗U_z;
%--------------------------------------------- Plot Result
figure (2); surf(H,B,U_MAU); shading interp; colormap(jet); colorbar; view(0,90); title('MAU′); xlabel('h'), ylabel('b')
%--------------------------------------------- Search for optimum
[Row,Col] = find(U_MAU = = max(max(U_MAU)));
B_H_O = [b_l+(Row-1)∗interval h_l+(Col-1)∗interval]
W_Z_O = [rou∗B_H_O(1)∗B_H_O(2)∗l∗9.8 P∗l ˆ 3/(3∗E∗(B_H_O(1)∗B_H_O(2) ˆ 3/12))]
%end
Script 4: solving the constrained problem of the cantilever beam example in Section 2.6.1.2.
%Unit: mm,N,MPa
clear all
%------------- Parameters that can be adjusted
k_w=0.5; k_z=0.5;           %%%%% Weighting between attributes
r_w = 0.0001;  r_z = −0.0001;    %%%%% Model risk; r > 0 risk adverse; r < 0 risk adverse; r∼0 risk neutral
slope = −0.01;            %%%%% Slope of the constraint utility function
k_C = 0.2;               %%%%% parameter for MAU with constraint
b_u = 60; b_l = 10;          %%%%% Upper and lower bounds of base
h_u = 80; h_l = 20;          %%%%% Upper and lower bounds of height
interval = 0.1;            %%%%% Solution search interval size
%------------- Define properties
l = 200;  P = 1000;    E = 69,000;  YS = 27.6;
rou = 2700∗10 ˆ (-9);  %Kg/mm3
%------------------------- An example; special case
b=60; h=80; I=b∗h ˆ 3/12; z=P∗l ˆ 3/(3∗E∗I); M=P∗l; stress = M∗(h/2)/I; w=rou∗b∗h∗l∗9.8;
%----------------------------------- Utility functions for individual attributes
w_l = rou∗b_u∗h_u∗l∗9.8;         w_u = rou∗b_l∗h_l∗l∗9.8;
z_l = P∗lˆ3/(3∗E∗(b_l∗h_l ˆ 3/12));     z_u = P∗l ˆ 3/(3∗E∗(b_u∗h_u ˆ 3/12));
stress_l = P∗l∗(h_u/2)/(b_u∗h_u ˆ 3/12);  stress_u = P∗l∗(h_l/2)/(b_l∗h_l ˆ 3/12);
 
x_w = w_l: 0.01:w_u;           x_z = z_l: 0.001:z_u;
s_w=(x_w-w_l)./(w_u-w_l);        s_z=(x_z-z_l)./(z_u-z_l);
u_w=(1-exp(-r_w.∗s_w))./(1-exp(-r_w));  u_z=(1-exp(-r_z.∗s_z))./(1-exp(-r_z)); 
 
figure (1);  % plot utility curves
subplot(1,3,1); plot(x_w,u_w); xlabel('Self weight'), ylabel('Utility')
subplot(1,3,2); plot(x_z,u_z); xlabel('Displacement'), ylabel('Utility')
x_stress = stress_l:0.1:stress_u;
u_stress = 1./(1 + exp((YS-x_stress)./slope));
subplot(1,3,3); plot(x_stress, u_stress);
xlabel('Stress'), ylabel('Utility'), axis([stress_l stress_u −0.1 1.1]);
%------------------------------------ Build MAU
K=(1-k_w-k_z)/(k_w∗k_z);
 
B = b_l:interval:b_u;
H = h_l:interval:h_u;
W = rou∗B.'∗H∗l∗9.8;
Z = P.∗l.ˆ 3./(3∗E∗(B.'∗H.ˆ 3/12));
 
S_w=(W-w_l)./(w_u-w_l);
S_z=(Z-z_l)./(z_u-z_l);
 
U_w=(1-exp(-r_w.∗S_w))./(1-exp(-r_w));
U_z=(1-exp(-r_w.∗S_z))./(1-exp(-r_w));
U_MAU = k_w∗U_w + k_z∗U_z + K∗k_w∗k_z∗U_w.∗U_z;  %MAU w/o constraint
%--------------------------------------------- Constraint
STRESS = 6.∗P.∗l./(B.'∗H.ˆ 2);
U_C = 1./(1 + exp((YS-STRESS)./slope));
U_MAU_C=(k_C+(1-k_C)∗U_MAU).∗U_C;   %MAU with constraint
%--------------------------------------------- Plot Result
figure (2); surf(H,B,STRESS); shading interp; colormap(jet); colorbar; view(0,90); title('Stress'), xlabel('h'), ylabel('b')
figure (3); surf(H,B,U_C); shadinginterp; colormap(jet); colorbar; view(0,90); title('Constraint'), xlabel('h'), ylabel('b')
figure (4); surf(H,B,U_MAU); shading interp; colormap(jet); colorbar; view(0,90); title('MAU′); xlabel('h'), ylabel('b')
figure (5); surf(H,B,U_MAU_C); shading interp; colormap(jet); colorbar; view(0,90); title('MAU with constraint'), xlabel('h'), ylabel('b')
%--------------------------------------------- Search for optimum
[Row,Col] = find(U_MAU = = max(max(U_MAU)));
[Row_C,Col_C] = find(U_MAU_C = = max(max(U_MAU_C)));
B_H_O = [b_l+(Row-1)∗interval h_l+(Col-1)∗interval]
W_Z_O = [rou∗B_H_O(1)∗B_H_O(2)∗l∗9.8 P∗l ˆ 3/(3∗E∗(B_H_O(1)∗B_H_O(2) ˆ 3/12))]
B_H_O_C = [b_l+(Row_C-1)∗interval h_l+(Col_C-1)∗interval]
W_Z_O_C = [rou∗B_H_O_C(1)∗B_H_O_C(2)∗l∗9.8 P∗l ˆ 3/(3∗E∗(B_H_O_C(1)∗B_H_O_C(2) ˆ 3/12))]
%end
Script 5: graphing the Pareto solution of Figure 2.27 in Section 2.6.2.3.
NN = 1,000,000;    % No. of points
%------------------
T_u = 6; T_l = 0.5;          %%%%% Upper and lower bounds of thickness
R_u = 38; R_l = 0;          %%%%% Upper and lower bounds of radius
%------------- Define properties
rou=0.283; L=100; St=32; P=4;
%------------------------
A = 0; B = 0;     % Count number of feasible and infeasible points
for ii = 1:NN
    TT = T_l + rand∗(T_u-T_l);
    RR = R_l + rand∗(R_u-R_l);
    WGT = rou∗pi∗((RR + TT) ˆ 2∗(L+2∗TT)-RR ˆ 2∗L);
    VOL = -pi∗RR ˆ 2∗L;
    stress = P∗RR/TT;
    if stress < St && 5∗TT ≤ RR && RR ≤ 40-TT
        A = A+1;
        Feasi(:,A) = [WGT,VOL];
    else
        B=B+1;
        Infeasi(:,B) = [WGT,VOL];
    end
end
figure (5)
hold on
plot(Feasi(1,:),Feasi(2,:),'b.','MarkerSize',8)    % Plot feasible solutions
%plot(Infeasi(1,:),Infeasi(2,:),'r.','MarkerSize',8)  % Plot infeasible solutions
%legend('Feasible','Infeasible')                             
xlabel('Self-weight (W)'), ylabel('Negative volume (-V)')
%-------------------------------------------- Plot Nash solution
T_Nash = 0.5:0.1:4.4;
R_Nash = 8.∗T_Nash;
WGT_Nash = rou.∗pi.∗((R_Nash + T_Nash).ˆ 2.∗(L+2.∗T_Nash)-R_Nash.ˆ 2.∗L);
VOL_Nash = -pi.∗R_Nash.ˆ 2.∗L;
plot(WGT_Nash, VOL_Nash,'k.','MarkerSize',20)
hold off
..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.
Reset