MrXuan_Bao2
Dinh Phuong
2025-10-17
#0 Add font, library
## Them font chu
library(showtext)## Loading required package: sysfonts## Loading required package: showtextdbfont_add("Times New Roman", "times.ttf") # Đường dẫn tới file font Times New Roman
font_add("Arial", "arial.ttf")
font_add("Times New Roman Bold", "timesbd.ttf")
showtext_auto()
# 1
library(ggplot2)
library(viridis)## Loading required package: viridisLitelibrary(ggthemes)
library(rsm)
# Bảng màu dùng chung
viridis_col_120 <- viridis(120, option = "D") # Mặc định: viridis (xanh tím → vàng)
viridis_col_full <- viridis(120, option = "turbo")
viridis_col_simple <- viridis_col_full[1:90] # Chọn từ vị trí 1 đến 90
my_col <- c("#0000FF", "#00FFFF", "#00FF00", "#FFFF00", "#FF8000", "#FF0000")#1 Chuẩn bị mô hình ##1.1. Mô hình QHTN Box Behnken
library(rsm)
# Box-Wilson Design bbd
QHTN <- bbd(
k = 3, # 3 yếu tố
n0 = 3, # 3 điểm trung tâm
block = FALSE, # No blocks
randomize = FALSE, # Not randomized
coding = list(
x1 ~ (KcTam_e - 25)/10,
x2 ~ (Gocluon_R - 5)/2,
x3 ~ (HsMasat_μ - 0.15)/0.05
)
)
print(QHTN)## run.order std.order KcTam_e Gocluon_R HsMasat_μ
## 1 1 1 15 3 0.15
## 2 2 2 35 3 0.15
## 3 3 3 15 7 0.15
## 4 4 4 35 7 0.15
## 5 5 5 15 5 0.10
## 6 6 6 35 5 0.10
## 7 7 7 15 5 0.20
## 8 8 8 35 5 0.20
## 9 9 9 25 3 0.10
## 10 10 10 25 7 0.10
## 11 11 11 25 3 0.20
## 12 12 12 25 7 0.20
## 13 13 13 25 5 0.15
## 14 14 14 25 5 0.15
## 15 15 15 25 5 0.15
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (KcTam_e - 25)/10
## x2 ~ (Gocluon_R - 5)/2
## x3 ~ (HsMasat_μ - 0.15)/0.05##1.2. Thí nghiệm thực tế và đưa dữ liệu vào R
## Ở đây mình chạy số liệu mô phỏng về LucF và Hệ số Ci
library(readxl)
data <- read_excel("~/1.Detai_Caohoc/MrXuan/Bao2/DataQHTN_ECAPPC_211025.xlsx", sheet = "BoxBehnken", col_types = c("numeric",
"numeric", "numeric", "numeric",
"numeric"))
head(data)## # A tibble: 6 × 5
## KCTam_e GocLuon_R HesoMS_f LucF_Y1 CiBD_Y2
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 15 3 0.15 124. 0.499
## 2 35 3 0.15 131. 0.340
## 3 15 7 0.15 116. 0.485
## 4 35 7 0.15 117. 0.578
## 5 15 5 0.1 111. 0.544
## 6 35 5 0.1 113. 0.457# Kiểm tra số dòng của QHTN và data_excel
if (nrow(data) != nrow(QHTN)) {
stop("Số dòng của dữ liệu Excel không khớp với số dòng của QHTN!")}
# Gắn cột y1 và y2 từ data vào QHTN
QHTN$LucF_Y1 <- data$LucF_Y1
QHTN$CiBD_Y2 <- data$CiBD_Y2
# In kết quả để kiểm tra
print(QHTN)## run.order std.order KcTam_e Gocluon_R HsMasat_μ LucF_Y1 CiBD_Y2
## 1 1 1 15 3 0.15 123.7659 0.4994021
## 2 2 2 35 3 0.15 130.7823 0.3400517
## 3 3 3 15 7 0.15 116.1473 0.4849841
## 4 4 4 35 7 0.15 116.7619 0.5778713
## 5 5 5 15 5 0.10 110.9146 0.5444828
## 6 6 6 35 5 0.10 113.0822 0.4565588
## 7 7 7 15 5 0.20 133.1021 0.6148683
## 8 8 8 35 5 0.20 133.7514 0.5768149
## 9 9 9 25 3 0.10 121.8491 0.3216295
## 10 10 10 25 7 0.10 107.1278 0.4771377
## 11 11 11 25 3 0.20 133.2195 0.3620474
## 12 12 12 25 7 0.20 130.8387 0.5030914
## 13 13 13 25 5 0.15 124.2060 0.4579854
## 14 14 14 25 5 0.15 125.6337 0.4834678
## 15 15 15 25 5 0.15 127.0610 0.5089501
##
## Data are stored in coded form using these coding formulas ...
## x1 ~ (KcTam_e - 25)/10
## x2 ~ (Gocluon_R - 5)/2
## x3 ~ (HsMasat_μ - 0.15)/0.05#2. Tính toán Hàm hồi quy ##2.1. Force F (ton)
#=============================================================================
# Box Behnken - Y1 Force F (ton)
#=============================================================================
ket_qua_model <- rsm((LucF_Y1) ~ FO(x1, x2, x3)+TWI(x1,x2,x3)+PQ(x1,x2,x3),
data = QHTN)
summary(ket_qua_model)##
## Call:
## rsm(formula = (LucF_Y1) ~ FO(x1, x2, x3) + TWI(x1, x2, x3) +
## PQ(x1, x2, x3), data = QHTN)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.63356 1.06778 117.6587 8.411e-10 ***
## x1 1.30600 0.65388 1.9973 0.1022901
## x2 -4.84264 0.65388 -7.4060 0.0007064 ***
## x3 9.74223 0.65388 14.8991 2.465e-05 ***
## x1:x2 -1.60043 0.92472 -1.7307 0.1440592
## x1:x3 -0.37958 0.92472 -0.4105 0.6984543
## x2:x3 3.08516 0.92472 3.3363 0.0206360 *
## x1^2 -2.15771 0.96248 -2.2418 0.0750438 .
## x2^2 -1.61149 0.96248 -1.6743 0.1549243
## x3^2 -0.76328 0.96248 -0.7930 0.4637111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9837, Adjusted R-squared: 0.9545
## F-statistic: 33.63 on 9 and 5 DF, p-value: 0.000602
##
## Analysis of Variance Table
##
## Response: (LucF_Y1)
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 960.54 320.18 93.6077 8.194e-05
## TWI(x1, x2, x3) 3 48.89 16.30 4.7649 0.06283
## PQ(x1, x2, x3) 3 25.86 8.62 2.5204 0.17199
## Residuals 5 17.10 3.42
## Lack of fit 3 13.03 4.34 2.1309 0.33523
## Pure error 2 4.08 2.04
##
## Stationary point of response surface:
## x1 x2 x3
## 1.735140 -3.565064 -1.254566
##
## Stationary point in original units:
## KcTam_e Gocluon_R HsMasat_μ
## 42.35139736 -2.13012879 0.08727168
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.5655325 -1.9867145 -3.1112960
##
## $vectors
## [,1] [,2] [,3]
## x1 0.2338790 0.8179920 -0.5255375
## x2 -0.6178802 -0.2922928 -0.7299240
## x3 -0.7506828 0.4954331 0.4370600# pred.r_squared
PRESS <- sum((residuals(ket_qua_model)/(1-lm.influence(ket_qua_model)$hat))^2)
SST <- sum((ket_qua_model$fitted.values + residuals(ket_qua_model) - mean(ket_qua_model$model[,1]))^2)
predR2 <- 1 - PRESS/SST
#cat("Predicted R-squared:", round(predR2, 4), "\n")
cat(sprintf("Std. Dev.: %.4f Mean: %.4f C.V.%%: %.2f%% Predicted R-squared: %.4f\n",
sigma(ket_qua_model), mean(ket_qua_model$model[,1]),
sigma(ket_qua_model)/mean(ket_qua_model$model[,1])*100, round(predR2, 4)))## Std. Dev.: 1.8494 Mean: 123.2162 C.V.%: 1.50% Predicted R-squared: 0.7932##2.2. Ci coefficient
#============================================================================
# Box Behnken - Y2 Ci coefficient
#============================================================================
ket_qua_modelY2_Full <- rsm((CiBD_Y2) ~ FO(x1, x2, x3)+TWI(x1,x2,x3)+PQ(x1,x2,x3),
data = QHTN)
summary(ket_qua_modelY2_Full)##
## Call:
## rsm(formula = (CiBD_Y2) ~ FO(x1, x2, x3) + TWI(x1, x2, x3) +
## PQ(x1, x2, x3), data = QHTN)
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4834678 0.0170112 28.4206 1.01e-06 ***
## x1 -0.0240551 0.0104172 -2.3092 0.068978 .
## x2 0.0649942 0.0104172 6.2391 0.001549 **
## x3 0.0321267 0.0104172 3.0840 0.027345 *
## x1:x2 0.0630594 0.0147321 4.2804 0.007860 **
## x1:x3 0.0124676 0.0147321 0.8463 0.436011
## x2:x3 -0.0036160 0.0147321 -0.2455 0.815863
## x1^2 0.0621571 0.0153337 4.0536 0.009790 **
## x2^2 -0.0700476 0.0153337 -4.5682 0.006012 **
## x3^2 0.0025563 0.0153337 0.1667 0.874129
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Multiple R-squared: 0.9577, Adjusted R-squared: 0.8817
## F-statistic: 12.59 on 9 and 5 DF, p-value: 0.006147
##
## Analysis of Variance Table
##
## Response: (CiBD_Y2)
## Df Sum Sq Mean Sq F value Pr(>F)
## FO(x1, x2, x3) 3 0.046680 0.0155601 17.9234 0.004168
## TWI(x1, x2, x3) 3 0.016580 0.0055267 6.3661 0.036856
## PQ(x1, x2, x3) 3 0.035099 0.0116997 13.4768 0.007857
## Residuals 5 0.004341 0.0008681
## Lack of fit 3 0.003042 0.0010140 1.5616 0.413323
## Pure error 2 0.001299 0.0006494
##
## Stationary point of response surface:
## x1 x2 x3
## 0.4472715 0.8402609 -6.7802361
##
## Stationary point in original units:
## KcTam_e Gocluon_R HsMasat_μ
## 29.4727146 6.6805218 -0.1890118
##
## Eigenanalysis:
## eigen() decomposition
## $values
## [1] 0.06977228 0.00219940 -0.07730586
##
## $vectors
## [,1] [,2] [,3]
## x1 0.97226548 -0.07353920 0.22201762
## x2 0.21815796 -0.05701057 -0.97424684
## x3 0.08430269 0.99566148 -0.03938625# pred.r_squared
PRESS <- sum((residuals(ket_qua_modelY2_Full)/(1-lm.influence(ket_qua_modelY2_Full)$hat))^2)
SST <- sum((ket_qua_modelY2_Full$fitted.values + residuals(ket_qua_modelY2_Full) - mean(ket_qua_modelY2_Full$model[,1]))^2)
predR2 <- 1 - PRESS/SST
#cat("Predicted R-squared:", round(predR2, 4), "\n")
cat(sprintf("Std. Dev.: %.4f Mean: %.4f C.V.%%: %.2f%% Predicted R-squared: %.4f\n",
sigma(ket_qua_modelY2_Full), mean(ket_qua_modelY2_Full$model[,1]),
sigma(ket_qua_modelY2_Full)/mean(ket_qua_modelY2_Full$model[,1])*100, round(predR2, 4)))## Std. Dev.: 0.0295 Mean: 0.4806 C.V.%: 6.13% Predicted R-squared: 0.4976### Fit Modified model ========================================================
ket_qua_modelY2 <- rsm(
CiBD_Y2 ~
x1 + x2 + x3 + # FO
x1:x2 + x1:x3 + # TWI (bỏ x2:x3)
I(x1^2) + I(x2^2), # PQ (bỏ x3^2)
data = QHTN
)## Warning in rsm(CiBD_Y2 ~ x1 + x2 + x3 + x1:x2 + x1:x3 + I(x1^2) + I(x2^2), : No FO() terms in model; cannot use RSM methods
## An 'lm' object has been returned.summary(ket_qua_modelY2)##
## Call:
## rsm(formula = CiBD_Y2 ~ x1 + x2 + x3 + x1:x2 + x1:x3 + I(x1^2) +
## I(x2^2), data = QHTN)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.027055 -0.013295 -0.001573 0.012339 0.029473
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.485041 0.012067 40.195 1.54e-09 ***
## x1 -0.024055 0.008881 -2.709 0.03026 *
## x2 0.064994 0.008881 7.318 0.00016 ***
## x3 0.032127 0.008881 3.617 0.00854 **
## I(x1^2) 0.061960 0.013034 4.754 0.00207 **
## I(x2^2) -0.070244 0.013034 -5.389 0.00102 **
## x1:x2 0.063059 0.012560 5.021 0.00153 **
## x1:x3 0.012468 0.012560 0.993 0.35395
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02512 on 7 degrees of freedom
## Multiple R-squared: 0.957, Adjusted R-squared: 0.914
## F-statistic: 22.25 on 7 and 7 DF, p-value: 0.0002823anova(ket_qua_modelY2)## Analysis of Variance Table
##
## Response: CiBD_Y2
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 0.004629 0.004629 7.3360 0.0302641 *
## x2 1 0.033794 0.033794 53.5547 0.0001602 ***
## x3 1 0.008257 0.008257 13.0851 0.0085404 **
## I(x1^2) 1 0.016748 0.016748 26.5410 0.0013214 **
## I(x2^2) 1 0.018327 0.018327 29.0438 0.0010202 **
## x1:x2 1 0.015906 0.015906 25.2068 0.0015294 **
## x1:x3 1 0.000622 0.000622 0.9853 0.3539501
## Residuals 7 0.004417 0.000631
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ====== THÊM PHẦN LACK OF FIT VÀ PURE ERROR ===================================
### Tính Lack of Fit và Pure Error cho mô hình modified
# 1. Lấy residual
resid <- residuals(ket_qua_modelY2)
# 2. Tính SSE (model error)
SSE <- sum(resid^2)
# 3. Tính Pure Error từ các điểm trùng (replicate) – dùng x1,x2,x3 thật
df_data <- QHTN
df_data$Y <- QHTN$CiBD_Y2
# Nhóm theo các điểm thiết kế trùng nhau
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionpure_error_data <- df_data %>%
group_by(x1, x2, x3) %>%
mutate(n = n()) %>%
filter(n > 1) %>% # chỉ tính các điểm lặp
summarise(SS = sum((Y - mean(Y))^2),
df = n() - 1)## `summarise()` has grouped output by 'x1', 'x2'. You can override using the
## `.groups` argument.SSPE <- sum(pure_error_data$SS)
df_PE <- sum(pure_error_data$df)
# 4. df tổng error của mô hình
n <- nrow(QHTN)
p <- length(coef(ket_qua_modelY2))
df_E <- n - p
# 5. Lack of fit
SSLOF <- SSE - SSPE
df_LOF <- df_E - df_PE
F_LOF <- (SSLOF / df_LOF) / (SSPE / df_PE)
p_LOF <- 1 - pf(F_LOF, df_LOF, df_PE)
cat("Lack of Fit vs Pure Error\n")## Lack of Fit vs Pure Errorcat(sprintf("SSLOF = %.4f df_LOF = %d\n", SSLOF, df_LOF))## SSLOF = 0.0031 df_LOF = 5cat(sprintf("SSPE = %.4f df_PE = %d\n", SSPE, df_PE))## SSPE = 0.0013 df_PE = 2cat(sprintf("F-value = %.4f p-value = %.6f\n\n", F_LOF, p_LOF))## F-value = 0.9605 p-value = 0.581217# ==============================================================================
# pred.r_squared
PRESS <- sum((residuals(ket_qua_modelY2)/(1-lm.influence(ket_qua_modelY2)$hat))^2)
SST <- sum((ket_qua_modelY2$fitted.values + residuals(ket_qua_modelY2) - mean(ket_qua_modelY2$model[,1]))^2)
predR2 <- 1 - PRESS/SST
#cat("Predicted R-squared:", round(predR2, 4), "\n")
cat(sprintf("Std. Dev.: %.4f Mean: %.4f C.V.%%: %.2f%% Predicted R-squared: %.4f\n",
sigma(ket_qua_modelY2), mean(ket_qua_modelY2$model[,1]),
sigma(ket_qua_modelY2)/mean(ket_qua_modelY2$model[,1])*100, round(predR2, 4)))## Std. Dev.: 0.0251 Mean: 0.4806 C.V.%: 5.23% Predicted R-squared: 0.7557#3. Đồ thị ##3.1. Đồ thị của Force F (Ton) ###3.1.1. Đồ thị Diagnostics
#11111111111111111111111111111111111111111111111111111111111111111111111111111
#Biểu đồ 1
#===============================
# NORMAL PROBABILITY PLOT – MÀU GIỐNG HỆT contour(rsm, image=TRUE)
#===============================
library(ggplot2)
library(showtext)
showtext_auto()
# 1. Externally Studentized Residuals
esr <- rstudent(ket_qua_model)
# 2. Tạo data frame + lấy giá trị Y1 để tô màu (hoặc Y2 nếu bạn muốn)
n <- length(esr)
df_plot <- data.frame(
ESR = sort(esr),
Theoretical= qnorm((1:n - 0.375)/(n + 0.25)), # Blom's formula
Y1 = QHTN$LucF_Y1[order(esr)] # Tô màu theo Lực F
# Y2 = QHTN$CiBD_Y2[order(esr)] # ← đổi thành Y2 nếu cần
)
# 3. Biểu đồ – màu giống hệt contour(rsm, image=TRUE)
F1 <- ggplot(df_plot, aes(x = ESR, y = Theoretical)) +
# Điểm được tô màu theo Y1, dùng đúng palette của rsm
geom_point(aes(fill = Y1),
shape = 22, size = 2, color = "black", stroke = 0.5) +
# Đường hồi quy đỏ kéo dài từ đầu đến cuối
geom_smooth(method = "lm", se = FALSE,
color = "red", linewidth = 1, fullrange = TRUE) +
# ←←←←←← QUAN TRỌNG: Dùng đúng palette terrain.colors của rsm
scale_fill_gradientn(
colours = viridis_col_simple, # ←←←←←← giống hệt contour(rsm, image=TRUE)
name = expression(bold("Lực F (Y1)"))
) +
# Trục x
scale_x_continuous(
name = "Externally Studentized Residuals",
limits = c(-3.2, 3.2), breaks = -3:3
) +
# Trục y: Normal % Probability chuẩn
scale_y_continuous(
name = "Normal % Probability",
breaks = qnorm(c(0.01,0.05,0.1,0.2,0.3,0.5,0.7,0.8,0.9,0.95,0.99)),
labels = c("1","5","10","20","30","50","70","80","90","95","99"),
limits = qnorm(c(0.005, 0.995))
) +
ggtitle("Normal Probability Plot of Residuals") +
theme_bw(base_family = "Times New Roman", base_size = 16) +
guides(fill = guide_colourbar(title.position = "top", title.hjust = 0.5))+
theme(
panel.grid.major.x = element_blank(), # Loại bỏ đường lưới dọc
panel.grid.minor.x = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.minor.y = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.major.y = element_blank(), # Loại bỏ đường lưới dọc
axis.ticks = element_line(linewidth = 0.5), #nét gạchight (tick marks) trên trục x và y
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5), # Thêm viền đen
axis.text = element_text(color = "black", size = 72), # Chữ trên trục màu đen
axis.title = element_text(color = "black", size = 72))+ # Tiêu đề trục màu đen
theme(plot.title = element_text(family = "Times New Roman Bold", size = 18L, hjust = 0.5),
axis.title.y = element_text(family = "Times New Roman", size = 36L, vjust = 2),
axis.title.x = element_text(family = "Times New Roman", size = 36L, , vjust = -2),
axis.text.y = element_text(family = "Times New Roman", size = 36L),
axis.text.x = element_text(family = "Times New Roman", size = 36L))+
theme(legend.title = element_text(family = "Times New Roman", size = 36L),
legend.text = element_text(family = "Times New Roman", size = 36L),
aspect.ratio = 1,
plot.margin = margin(10, 10, 10, 10))+
theme(legend.position = "none", #Xoá legend
plot.title = element_blank()) #Xoá title
# Hiển thị
print(F1)## `geom_smooth()` using formula = 'y ~ x'
#png("F1.png", width = 8, height = 8, units = "cm", res = 600)
#plot(F1, main = "My Base R Plot")
#dev.off()
#2222222222222222222222222222222222222222222222222222222222222222222222222222
library(ggplot2)
library(showtext)
showtext_auto()
model <- ket_qua_model
# ========== TỰ ĐỘNG TÍNH GIỚI HẠN ĐÚNG NHƯ DESIGN-EXPERT ==========
n <- nrow(QHTN)
p <- length(coef(model))
df <- n - p - 1
alpha <- 0.05
critical_t <- qt(1 - alpha/(2*n), df) # Bonferroni critical value
lower_limit <- -critical_t
upper_limit <- +critical_t
cat(sprintf("Số quan sát n = %d, số tham số p = %d → df = %d\n", n, p, df))## Số quan sát n = 15, số tham số p = 10 → df = 4cat(sprintf("Giới hạn Bonferroni (α=0.05): ± %.4f\n", critical_t))## Giới hạn Bonferroni (α=0.05): ± 6.2541# Ví dụ: bạn sẽ thấy ±5.8, ±6.1, hoặc ±6.25 tùy model → khớp 100% với Design-Expert!
# =================================================================
esr <- rstudent(model)
run_order <- 1:length(esr)
df_run <- data.frame(
Run = run_order,
ESR = esr,
Y1 = QHTN$LucF_Y1
)
F2 <- ggplot(df_run, aes(x = Run, y = ESR)) +
geom_hline(yintercept = 0, color = "black", linewidth = 0.8) +
# Đường giới hạn TỰ ĐỘNG, CHÍNH XÁC NHƯ DESIGN-EXPERT
geom_hline(yintercept = lower_limit, color = "red", linetype = "solid", linewidth = 1) +
geom_hline(yintercept = upper_limit, color = "red", linetype = "solid", linewidth = 1) +
# Ghi chú giá trị giới hạn lên đồ thị (rất chuyên nghiệp)
annotate("text", x = 1, y = upper_limit,
label = sprintf("%.3f", upper_limit),
vjust = -0.8, hjust = 0, color = "red", size = 12, family = "Times New Roman") +
annotate("text", x = 1, y = lower_limit,
label = sprintf("%.3f", lower_limit),
vjust = 1.8, hjust = 0, color = "red", size = 12, family = "Times New Roman") +
geom_point(aes(fill = Y1), shape = 22, size = 2, color = "black", stroke = 0.5) +
geom_line(color = "black", linewidth = 0.5) +
scale_fill_gradientn(colours = viridis_col_simple,
name = expression(bold("Lực F (Y1)"))) +
scale_x_continuous(name = "Run Number",
breaks = seq(1, max(run_order), by = 2), # chỉ hiện số lẻ: 1,3,5,...,15
minor_breaks = 1:max(run_order), # vẫn có gạch ngắn ở tất cả các số
expand = c(0.02, 0))+ # sát viền đẹp hơn
scale_y_continuous(name = "Externally Studentized Residuals",
limits = c(-8.2, 8.2),
#labels = function(x) sprintf("%.2f", x),
breaks = round(seq(-8, 8, 2), 1)) +
ggtitle("Residuals vs. Run Order\n(Externally Studentized Residuals with Bonferroni Limits)") +
theme_bw(base_family = "Times New Roman", base_size = 16) +
theme(
panel.grid.major.x = element_blank(), # Loại bỏ đường lưới dọc
panel.grid.minor.x = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.minor.y = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.major.y = element_blank(), # Loại bỏ đường lưới dọc
axis.ticks = element_line(linewidth = 0.5), #nét gạchight (tick marks) trên trục x và y
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5), # Thêm viền đen
axis.text = element_text(color = "black", size = 12), # Chữ trên trục màu đen
axis.title = element_text(color = "black", size = 12))+ # Tiêu đề trục màu đen
theme(plot.title = element_text(family = "Times New Roman Bold", size = 18L, hjust = 0.5),
axis.title.y = element_text(family = "Times New Roman", size = 36L, vjust = 2),
axis.title.x = element_text(family = "Times New Roman", size = 36L, , vjust = -2),
axis.text.y = element_text(family = "Times New Roman", size = 36L),
axis.text.x = element_text(family = "Times New Roman", size = 36L))+
theme(legend.title = element_text(family = "Times New Roman", size = 36L),
legend.text = element_text(family = "Times New Roman", size = 36L),
aspect.ratio = 1,
plot.margin = margin(10, 10, 10, 10))+
theme(legend.position = "none", #Xoá legend
plot.title = element_blank()) #Xoá title
print(F2)
#png("F2.png", width = 8, height = 8, units = "cm", res = 600)
#plot(F2, main = "My Base R Plot")
#dev.off()
#3333333333333333333333333333333333333333333333333333333333333333333333333333
# Biểu đồ 3
# PREDICTED vs ACTUAL
# ===============================
library(ggplot2)
library(showtext)
showtext_auto()
model <- ket_qua_model
actual <- QHTN$LucF_Y1 # hoặc CiBD_Y2 nếu bạn đang làm Y2
predicted <- predict(model)
df_pred <- data.frame(
Actual = actual,
Predicted = predicted,
Y1 = actual # tô màu theo giá trị thực tế
)
# ---- Tự động xác định min–max và tạo breaks đẹp ----
range_val <- range(c(df_pred$Actual, df_pred$Predicted)) # min và max chung
min_val <- floor(range_val[1] / 5) * 5 # làm tròn xuống
max_val <- ceiling(range_val[2] / 5) * 5 # làm tròn lên
break_seq <- seq(min_val, max_val, by = 5) # bạn thay 5 → 10, 20… tùy thích
# ---- Biểu đồ hoàn hảo ----
F3 <- ggplot(df_pred, aes(x = Actual, y = Predicted)) +
# Đường hồi quy + vùng tin cậy 95%
geom_abline(intercept = 0, slope = 1, color = "black", linewidth = 0.5, linetype = "solid") +
# Các điểm tô màu giống contour(rsm)
geom_point(aes(fill = Y1), shape = 22, size = 2, color = "black", stroke = 0.5) +
# Palette terrain.colors giống hệt 2 đồ thị trước
scale_fill_gradientn(colours = viridis_col_simple,
name = expression(bold("Force F (ton)"))) +
# Trục x, y: tự động min–max, nhưng breaks tùy chỉnh
scale_x_continuous(name = "Actual",
limits = c(min_val, max_val),
breaks = break_seq) +
scale_y_continuous(name = "Predicted",
limits = c(min_val, max_val),
breaks = break_seq) +
ggtitle("Predicted vs. Actual Values") +
theme_bw(base_family = "Times New Roman", base_size = 12) +
theme(
panel.grid.major.x = element_blank(), # Loại bỏ đường lưới dọc
panel.grid.minor.x = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.minor.y = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.major.y = element_blank(), # Loại bỏ đường lưới dọc
axis.ticks = element_line(linewidth = 0.5), #nét gạchight (tick marks) trên trục x và y
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5), # Thêm viền đen
axis.text = element_text(color = "black", size = 12), # Chữ trên trục màu đen
axis.title = element_text(color = "black", size = 12))+ # Tiêu đề trục màu đen
theme(plot.title = element_text(family = "Times New Roman Bold", size = 12L, hjust = 0.5),
axis.title.y = element_text(family = "Times New Roman", size = 36L, vjust = 2),
axis.title.x = element_text(family = "Times New Roman", size = 36L, , vjust = -2),
axis.text.y = element_text(family = "Times New Roman", size = 36L),
axis.text.x = element_text(family = "Times New Roman", size = 36L),
aspect.ratio = 1,
plot.margin = margin(10, 10, 10, 10))+
theme(legend.position = "right", # hoặc "bottom" nếu muốn nằm dưới
legend.title = element_text(family = "Times New Roman", size = 36L, hjust = 0.5), # tiêu đề legend
legend.text = element_text(family = "Times New Roman", size = 36L), # số trên thanh màu
legend.key.height = unit(0.5, "cm"), # chiều cao thanh màu (tăng/giảm thoải mái)
legend.key.width = unit(0.5, "cm"), # độ dày thanh màu
#legend.background = element_rect(fill = "white", color = "black", linewidth = 0.5), # viền legend
legend.margin = margin(10, 10, 10, 10))+ # khoảng cách bên trong legend
theme(plot.title = element_blank()) #Xoá title
# Xem trước
print(F3)
#png("F3.png", width = 16, height = 8, units = "cm", res = 600)
#plot(F3, main = "My Base R Plot")
#dev.off()###3.1.2. Đồ thị One Factor
#44444444444444444444444444444444444444444444444444444444444444444444444444444
# Biểu đồ 4
# ===============================
# ONE FACTOR PLOT: Ảnh hưởng của R → Lực F
# Trục Y tự động min–max + tùy chỉnh khoảng cách chia (ví dụ: mỗi 10)
# ===============================
library(ggplot2)
library(dplyr)
library(showtext)
showtext_auto()
model <- ket_qua_model
# ---- 1. Tạo lưới dự báo cho R từ 3 → 7 mm (giữ x1 = x3 = 0) ----
R_real <- seq(3, 7, length.out = 200)
grid <- data.frame(
x1 = 0,
x2 = (R_real - 5)/2, # chuyển sang coded
x3 = 0
)
# Dự báo + 95% Confidence Interval
pred <- suppressWarnings(
predict(model, newdata = grid, interval = "confidence", level = 0.95)
)
grid_df <- data.frame(
R_real = R_real,
fit = pred[,"fit"],
lwr = pred[,"lwr"],
upr = pred[,"upr"]
)
# ---- 2. Lấy đúng 3 điểm trung tâm thực nghiệm (R = 5 mm) ----
center_points <- data %>%
filter(GocLuon_R == 5 & KCTam_e == 25 ) %>%
select(GocLuon_R, LucF_Y1)
# ---- 3. TỰ ĐỘNG xác định min–max trục Y + tùy chỉnh breaks ----
all_y_values <- c(grid_df$lwr, grid_df$upr, center_points$LucF_Y1)
min_y <- floor(min(all_y_values) / 10) * 10 # làm tròn xuống đến bội số 10
max_y <- ceiling(max(all_y_values) / 10) * 10 # làm tròn lên đến bội số 10
# TỰ DO CHỈNH khoảng cách chia tại đây:
break_step <- 10 # ← Thay 10 → 20, 50, 5… tùy bạn
y_breaks <- seq((min_y)-10, max_y, by = break_step)
# ---- 4. Biểu đồ đẹp hơn cả Design-Expert ----
F4 <- ggplot() +
# Vùng 95% Confidence Interval
geom_ribbon(data = grid_df,
aes(x = R_real, ymin = lwr, ymax = upr),
fill = "gray", alpha = 0.3) +
# Đường dự báo trung bình
geom_line(data = grid_df,
aes(x = R_real, y = fit),
color = "black", linewidth = 0.5) +
# Chỉ 3 điểm đỏ ở tâm (R = 5)
geom_point(data = center_points,
aes(x = GocLuon_R, y = LucF_Y1),
color = "darkred", size = 1.2, shape = 19) +
# Trục X
scale_x_continuous(name = "Corner radius R (mm)",
breaks = seq(3, 7, by = 1)) +
# Trục Y: tự động min–max + bạn tùy chỉnh khoảng cách chia
scale_y_continuous(name = "Force F (ton)",
limits = c((min_y)-10, max_y),
breaks = y_breaks) + # ← Tùy chỉnh ở đây!
ggtitle("One Factor Plot: Effect of R on Lực F\n(x1 = 0, x3 = 0)") +
# Theme đồng bộ tuyệt đối với tất cả biểu đồ trước
theme_bw(base_family = "Times New Roman", base_size = 12) +
theme(
panel.grid.major.x = element_blank(), # Loại bỏ đường lưới dọc
panel.grid.minor.x = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.minor.y = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.major.y = element_blank(), # Loại bỏ đường lưới dọc
axis.ticks = element_line(linewidth = 0.5), #nét gạchight (tick marks) trên trục x và y
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5), # Thêm viền đen
axis.text = element_text(color = "black", size = 12), # Chữ trên trục màu đen
axis.title = element_text(color = "black", size = 12))+ # Tiêu đề trục màu đen
theme(plot.title = element_text(family = "Times New Roman Bold", size = 18L, hjust = 0.5),
axis.title.y = element_text(family = "Times New Roman", size = 36L, vjust = 2),
axis.title.x = element_text(family = "Times New Roman", size = 36L, , vjust = -2),
axis.text.y = element_text(family = "Times New Roman", size = 36L),
axis.text.x = element_text(family = "Times New Roman", size = 36L))+
theme(legend.title = element_text(family = "Times New Roman", size = 36L),
legend.text = element_text(family = "Times New Roman", size = 36L),
aspect.ratio = 1,
plot.margin = margin(10, 10, 10, 10))+
theme(legend.position = "none", #Xoá legend
plot.title = element_blank()) #Xoá title
# Xem trước
print(F4)
#png("F4.png", width = 8, height = 8, units = "cm", res = 600)
#plot(F4, main = "My Base R Plot")
#dev.off()
#5555555555555555555555555555555555555555555555555555555555555555555555555555
# Biểu đồ 5
# ===============================
# Biểu đồ One Factor: Ảnh hưởng của Hệ số ma sát μ → Lực F
# Hoàn toàn giống format F4 của bạn (8x8 cm, font 72, không title/legend)
# ===============================
library(ggplot2)
library(dplyr)
library(showtext)
showtext_auto()
model <- ket_qua_model # hoặc ket_qua_model_opt
# ---- 1. Tạo lưới dự báo cho μ từ 0.10 → 0.20 (giữ x1=x2=0) ----
mu_real <- seq(0.10, 0.20, length.out = 200)
grid <- data.frame(
x1 = 0,
x2 = 0,
x3 = (mu_real - 0.15)/0.05 # chuyển sang coded: (μ - 0.15)/0.05
)
# Dự báo + 95% Confidence Interval
pred <- suppressWarnings(
predict(model, newdata = grid, interval = "confidence", level = 0.95)
)
grid_df <- data.frame(
mu_real = mu_real,
fit = pred[,"fit"],
lwr = pred[,"lwr"],
upr = pred[,"upr"]
)
# 3. LẤY ĐIỂM TRUNG TÂM – DÙNG round() ĐỂ TRÁNH LỖI LÀM TRÒN
center_points <- data %>%
mutate(mu_round = round(HesoMS_f, 3)) %>% # làm tròn 3 chữ số
filter(mu_round == 0.150 & KCTam_e == 25) %>% # chắc chắn bắt được 0.15
select(mu_real = HesoMS_f, LucF_Y1) # đổi tên cột cho khớp
# Kiểm tra nhanh: phải ra đúng 3 dòng
print(center_points)## # A tibble: 3 × 2
## mu_real LucF_Y1
## <dbl> <dbl>
## 1 0.15 124.
## 2 0.15 126.
## 3 0.15 127.# ---- 3. TỰ ĐỘNG min–max trục Y + chia mỗi 10 ----
all_y <- c(grid_df$lwr, grid_df$upr, center_points$LucF_Y1)
min_y <- floor(min(all_y) / 10) * 10
max_y <- ceiling(max(all_y) / 10) * 10
break_step <- 10
y_breaks <- seq(min_y - 10, max_y, by = break_step)
# ---- 4. Biểu đồ giống hệt F4 của bạn ----
F5 <- ggplot() +
geom_ribbon(data = grid_df,
aes(x = mu_real, ymin = lwr, ymax = upr),
fill = "gray", alpha = 0.3) +
geom_line(data = grid_df,
aes(x = mu_real, y = fit),
color = "black", linewidth = 0.5) +
geom_point(data = center_points,
aes(x = mu_real, y = LucF_Y1),
color = "darkred", size = 1.2, shape = 19) +
scale_x_continuous(name = "Friction coefficient μ",
breaks = seq(0.10, 0.20, by = 0.02)) +
scale_y_continuous(name = "Force F (ton)",
limits = c(min_y - 10, max_y),
breaks = y_breaks) +
theme_bw(base_family = "Times New Roman", base_size = 12) +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5),
axis.ticks = element_line(linewidth = 0.5),
axis.text = element_text(color = "black", size = 36L, family = "Times New Roman"),
axis.title = element_text(color = "black", size = 36L, family = "Times New Roman"),
axis.title.y = element_text(vjust = 2),
axis.title.x = element_text(vjust = -2),
plot.margin = margin(10, 10, 10, 10),
legend.position = "none",
plot.title = element_blank()
)
# Xem trước
print(F5)
# Lưu file đúng chuẩn bạn đang dùng
#png("F5.png", width = 8, height = 8, units = "cm", res = 600)
#print(F5)
#dev.off()###3.1.3. Đồ thị Contour Plot
#666666666666666666666666666666666666666666666666666666666666666666666666666
#png("F6.png", width = 8, height = 8, units = "cm", res = 600)
library(viridis)
# Tạo dải 120 màu (giống col_dx_120 của bạn)
viridis_col_120 <- viridis(120, option = "D") # Mặc định: viridis "D" (xanh tím → vàng) "turbo" là tím-đỏ
cividis_col_120 <- viridis(120, option = "E") # Xanh → vàng (colorblind-safe nhất)
my_col <- c("#0000FF", "#00FFFF", "#00FF00", "#FFFF00", "#FF8000", "#FF0000")
col_dx_120 <- colorRampPalette(my_col)(120)
viridis_col_full <- viridis(120, option = "turbo")
viridis_col_simple <- viridis_col_full[1:90] # Chọn từ vị trí 1 đến 90
par(family = "Times New Roman", ps = 36, cex = 1,
mar = c(4, 3.5, 1, 0.5), oma = c(0, 0, 0, 0), mgp = c(3, 0.6, 0), las = 1, tcl = -0.3) # Muốn export ra 600Dpi thì ps=72
# Tạo lưới đều
R_vals <- seq(3, 7, length.out = 300)
mu_vals <- seq(0.1, 0.2, length.out = 300)
# Mở rộng thêm 0.15 ở mỗi đầu → tạo đoạn thừa rõ rệt
xlim <- c(2.9, 7.1) # thay vì c(3,7)
ylim <- c(0.097, 0.203) # thay vì c(0.10, 0.20)
grid <- expand.grid(x2 = (R_vals - 5)/2, x3 = (mu_vals - 0.15)/0.05)
grid$x1 <- 0
grid$yhat <- predict(ket_qua_model, grid)
zmat <- matrix(grid$yhat, 300, 300)
# Vẽ nền gradient đúng màu bạn muốn
image(R_vals, mu_vals, zmat, col = viridis_col_simple, xlim = xlim, ylim = ylim,
xlab = "", ylab = "", axes = TRUE, useRaster = TRUE)
# Vẽ đường đồng mức + nhãn số
contour(R_vals, mu_vals, zmat, add = TRUE,
col = "black", lwd = 1, labcex = 1, drawlabels = TRUE, method = "flattest")
# Trục + nhãn
# Vẽ trục (chỉ vạch, không chữ)
#axis(1, tcl = -0.3, mgp = c(3, 0.6, 0))
#axis(2, tcl = -0.3, mgp = c(3, 0.8, 0), las = 1)
# THÊM NHÃN TRỤC ĐÚNG NHƯ BẠN MUỐN (đẹp, điều chỉnh được line, cex, font)
title(xlab = "Corner radius R (mm)",
line = 1.6, # chỉnh độ cao trục X ở đây
cex.lab = 1, # đúng 12 pt
family = "Times New Roman")
title(ylab = "Friction coefficient μ",
line = 2.4, # chỉnh độ cao trục Y ở đây
cex.lab = 1, # đúng 12 pt
family = "Times New Roman")
# Tiêu đề slice (bạn tự đặt, không còn "Slice at x1 = 0")
mtext("Slice at Distance e = 25 mm",
family = "Times New Roman",
side = 1, line = 2.8, cex = 1, font = 1) # font = 2 → in đậm
#dev.off()###3.1.4. Đồ thị 3D Surface Plot
#===================================================================CHUẨN MẪU
#777777777777777777777777777777777777777777777777777777v22222222222222222222
## Biểu đồ 3D Surface ảnh hưởng e & R tới Ci
#png("F7v2.png", width = 8, height = 8, units = "cm", res = 600)
viridis_col_full <- viridis(120, option = "turbo")
viridis_col_simple <- viridis_col_full[1:90] # Chọn từ vị trí 1 đến 90
par(
family = "Times New Roman",
ps =36,
cex = 1,
mar = c(2.6, 2.6, 0.1, 0.1), #(Dưới, Trái, Trên, Phải)
mgp = c(0.5, 0.02, 0.7) #(title, tick lable, tick mark)
)
# Tạo lưới đều
e_vals <- seq(15, 35, length.out = 300)
R_vals <- seq(3, 7, length.out = 300)
mu_vals <- seq(0.1, 0.2, length.out = 300)
grid <- expand.grid(x3 = (mu_vals - 0.15)/0.05, x2 = (R_vals - 5)/2)
grid$x1 <- 0
grid$yhat <- predict(ket_qua_model, grid)
zmat <- matrix(grid$yhat, 300, 300)
min_z <- min(zmat)
max_z <- max(zmat)
# Tạo màu trống phía trên, dưới
Color_Tren <- ((140-max_z)/(max_z - min_z))*90
Color_Duoi <-((min_z-90)/(max_z - min_z))*90
# Thêm các màu trống cả phía trước và phía sau
viridis_col_extended <- c(
rep(viridis_col_simple[1], Color_Duoi), # Thêm 36 màu giống màu đầu tiên (phía dưới)
viridis_col_simple, # Dải màu chính từ thấp đến cao
rep(viridis_col_simple[90], Color_Tren) # Thêm 36 màu giống màu cuối cùng (phía trên)
)
# Tạo các mức contour cách nhau 0.05
levels_05 <- seq(
floor(min_z/0.05)*0.05,
ceiling(max_z/0.05)*0.05,
by = 0.05)
# Vẽ đồ thị không có trục
res <- persp(
ket_qua_model,
~ x2 + x3,
at = c(x1 = 0),
zlim = c(90, 140), # ← thêm khoảng trống !!!
axes = FALSE,
#sub = NULL, # ẩn Slice at
col = viridis_col_extended,
contours = "colors",
#contours = list(col = "colors", levels = my_levels), # <<< QUAN TRỌNG
color.legend = TRUE,
expand = 1,
theta = -28,
phi = 20
)## Warning in persp.default(dat$x, dat$y, dat$z, zlim = dat$zlim, theta = theta, :
## "color.legend" is not a graphical parameter## Warning in persp.default(dat$x, dat$y, dat$z, xlab = dat$labs[1], ylab =
## dat$labs[2], : "color.legend" is not a graphical parameter## Warning in title(sub = dat$labs[5], ...): "axes" is not a graphical parameter## Warning in title(sub = dat$labs[5], ...): "color.legend" is not a graphical
## parameter## Warning in title(sub = dat$labs[5], ...): "expand" is not a graphical parameter# Tạo lưới và các giá trị tick
x2_seq_coded <- seq(min(QHTN$x2), max(QHTN$x2), length.out = 80)
x3_seq_coded <- seq(min(QHTN$x3), max(QHTN$x3), length.out = 80)
#e_seq_coded <- seq(min(data$KCTam_e), max(data$KCTam_e), length.out = 80)
R_seq_coded <- seq(min(data$GocLuon_R), max(data$GocLuon_R), length.out = 80)
f_seq_coded <- seq(min(data$HesoMS_f), max(data$HesoMS_f), length.out = 80)
grid <- expand.grid(x2 = x2_seq_coded, x3 = x3_seq_coded, x1 = 0)
grid_real <- expand.grid(x1 = 25, x2 = R_seq_coded, x3 = f_seq_coded)
grid$F_pred <- predict(ket_qua_model, newdata = grid)
z_mat <- matrix(grid$F_pred, nrow = length(x2_seq_coded), ncol = length(x3_seq_coded))
x_ticks1 <- c("–","–","–","–","–")
x_ticks2 <- c("3","4","5","6","7")
y_ticks1 <- c("–","–","–","–","–","–")
y_ticks2 <- c("0.10","0.12","0.14","0.16","0.18","0.20")
z_ticks1 <- c("","–","–","–","–","–")
z_ticks2 <- c("","100","110","120","130","140")
par(xpd = NA)
#################################################### TRỤC X
###X11 Vẽ tick marks cho trục x2 với phép xoay toàn bộ trục
rotation_angle <- -73 # Góc xoay của toàn bộ trục x2 (độ)
angle_rad <- rotation_angle * pi / 180
offset1 <- 0.403
for(i in seq_along(x_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base11_x <- 0.247 #giảm là lên
base11_y <- -0.239 #giảm là sang trái
tick_offset <- (i/length(x_ticks1)) * offset1
x_before <- base11_x
y_before <- base11_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = x_ticks1[i],
cex = 1,
srt = rotation_angle)
}
###X12 Vẽ tick labels cho trục x2 với phép xoay toàn bộ trục
for(i in seq_along(x_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base12_x <- base11_x +0.03 #giảm là lên
base12_y <- base11_y #giảm là sang trái
tick_offset <- (i/length(x_ticks1)) * offset1
x_before <- base12_x
y_before <- base12_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = x_ticks2[i],
cex = 1,
srt = 0)
}
############ CÁCH 2: VẼ GẠCH TICK MARK BẰNG ĐIỂM
#for(i in seq_along(x_ticksmark)) {
# Tọa độ gốc của tick trước khi xoay
#tick_offset <- (i/length(x_ticksmark)) * 0.4
#x_before <- base_x
#y_before <- base_y + tick_offset
# Tính tọa độ sau khi xoay
#rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
#rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
#Vẽ tick mark tại vị trí đã xoay
#tmx21=rotated_x-0.062 #giảm là sang trái
#tmx22=rotated_y+0.01 #giảm là xuống
#segments(tmx21, tmx22,
#tmx21+0.008, tmx22-0.012,
#col = "black", lwd = 1.2)}
#################################################### TRỤC Y
###Y21 Vẽ tick marks cho trục x3 với phép xoay toàn bộ trục
rotation_angle <- 37 # Góc xoay của toàn bộ trục x3 (độ)
angle_rad <- rotation_angle * pi / 180
offset2 <- 0.236
for(i in seq_along(y_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base21_x <- -0.239 #giảm là sang trái
base21_y <- -0.2115 #giảm là xuống
tick_offset <- (i/length(y_ticks1)) * offset2
x_before <- base21_x
y_before <- base21_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = y_ticks1[i],
cex = 1,
srt = rotation_angle)
}
###Y22 Vẽ tick labels cho trục x3 với phép xoay toàn bộ trục
for(i in seq_along(y_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base22_x <- base21_x -0.04 #giảm là sang trái
base22_y <- base21_y +0.01 #giảm là xuống
tick_offset <- (i/length(y_ticks1)) * offset2
x_before <- base22_x
y_before <- base22_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay CÁCH 2
#segments(rotated_x - 0.05, rotated_y,
#rotated_x + 0.05, rotated_y,
#col = "black", lwd = 1.2)
# Vẽ tick label tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = y_ticks2[i],
cex = 1,
srt = rotation_angle-37)
}
#################################################### TRỤC Z
###Z31 Vẽ tick marks cho trục z với phép xoay toàn bộ trục
rotation_angle <- 4.2 # Góc xoay của toàn bộ trục z (độ)
angle_rad <- rotation_angle * pi / 180
offset3 <- 0.3568
for(i in seq_along(z_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base31_x <- -0.217 #giảm là sang trái
base31_y <- -0.163 #giảm là xuống
tick_offset <- (i/length(z_ticks1)) * offset3
x_before <- base31_x
y_before <- base31_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = z_ticks1[i],
cex = 1,
srt = rotation_angle)
}
###Z32 Vẽ tick label cho trục z với phép xoay toàn bộ trục
for(i in seq_along(z_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base32_x <- base31_x -0.04 #giảm là sang trái
base32_y <- base31_y #giảm là xuống
tick_offset <- (i/length(z_ticks1)) * offset3
x_before <- base32_x
y_before <- base32_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay CÁCH 2
#segments(rotated_x - 0.05, rotated_y,
#rotated_x + 0.05, rotated_y,
#col = "black", lwd = 1.2)
# Vẽ tick label tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = z_ticks2[i],
cex = 1,
srt = 0)
}
# Thêm nhãn cho các trục
par(xpd = NA, srt = -55)
text(-0.25, -0.26, "Friction coefficient μ")
par(xpd = NA, srt = 20)
text(0.13, -0.3, "Corner radius R (mm)")
par(xpd = NA, srt = -85)
text(-0.33, 0.02, "Force F (ton)")
#dev.off()##3.2. Đồ thị của Ci coefficient ###3.2.1. Đồ thị Diagnostics
#111111111111111111111111111111111111111111111111111111111111111111111111111
#Biểu đồ 1
#===============================
# NORMAL PROBABILITY PLOT – MÀU GIỐNG HỆT contour(rsm, image=TRUE)
#===============================
library(ggplot2)
library(showtext)
showtext_auto()
# 1. Externally Studentized Residuals
esr <- rstudent(ket_qua_modelY2)
# 2. Tạo data frame + lấy giá trị Y1 để tô màu (hoặc Y2 nếu bạn muốn)
n <- length(esr)
df_plot <- data.frame(
ESR = sort(esr),
Theoretical= qnorm((1:n - 0.375)/(n + 0.25)), # Blom's formula
#Y1 = QHTN$LucF_Y1[order(esr)] # Tô màu theo Lực F
Y2 = QHTN$CiBD_Y2[order(esr)] # ← đổi thành Y2 nếu cần
)
# 3. Biểu đồ – màu giống hệt contour(rsm, image=TRUE)
F21 <- ggplot(df_plot, aes(x = ESR, y = Theoretical)) +
# Điểm được tô màu theo Y2, dùng đúng palette của rsm
geom_point(aes(fill = Y2),
shape = 22, size = 2, color = "black", stroke = 0.5) +
# Đường hồi quy đỏ kéo dài từ đầu đến cuối NORMAL PROBABILITY PLOT là đường best fit
geom_smooth(method = "lm", se = FALSE,
color = "red", linewidth = 1, fullrange = TRUE) +
#geom_abline(intercept = 0, slope = 1, color = "black", linewidth = 0.5, linetype = "solid") +
# Các điểm tô màu giống contour(rsm)
# ←←←←←← QUAN TRỌNG: Dùng đúng palette terrain.colors của rsm
scale_fill_gradientn(
colours = viridis_col_simple, # ←←←←←← giống hệt contour(rsm, image=TRUE)
name = expression(bold("Ci (Y2)"))
) +
# Trục x được chia tự động với định dạng số thập phân 2 chữ số
#scale_x_continuous(
#name = "Externally Studentized Residuals",
#labels = function(x) sprintf("%.2f", x),
#n.breaks = 6 # Có thể thêm để gợi ý số lượng chia, nhưng không bắt buộc) +
# Trục x
scale_x_continuous(
name = "Externally Studentized Residuals",
#labels = function(x) sprintf("%.2f", x),
limits = c(-2.2, 3.2), breaks = -2:3) +
# Trục y: Normal % Probability chuẩn
scale_y_continuous(
name = "Normal % Probability",
breaks = qnorm(c(0.01,0.05,0.1,0.2,0.3,0.5,0.7,0.8,0.9,0.95,0.99)),
labels = c("1","5","10","20","30","50","70","80","90","95","99"),
limits = qnorm(c(0.005, 0.995))
) +
ggtitle("Normal Probability Plot of Residuals") +
theme_bw(base_family = "Times New Roman", base_size = 16) +
guides(fill = guide_colourbar(title.position = "top", title.hjust = 0.5))+
theme(
panel.grid.major.x = element_blank(), # Loại bỏ đường lưới dọc
panel.grid.minor.x = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.minor.y = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.major.y = element_blank(), # Loại bỏ đường lưới dọc
axis.ticks = element_line(linewidth = 0.5), #nét gạchight (tick marks) trên trục x và y
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5), # Thêm viền đen
axis.text = element_text(color = "black", size = 72), # Chữ trên trục màu đen
axis.title = element_text(color = "black", size = 72))+ # Tiêu đề trục màu đen
theme(plot.title = element_text(family = "Times New Roman Bold", size = 18L, hjust = 0.5),
axis.title.y = element_text(family = "Times New Roman", size = 36L, vjust = 2),
axis.title.x = element_text(family = "Times New Roman", size = 36L, , vjust = -2),
axis.text.y = element_text(family = "Times New Roman", size = 36L),
axis.text.x = element_text(family = "Times New Roman", size = 36L))+
theme(legend.title = element_text(family = "Times New Roman", size = 36L),
legend.text = element_text(family = "Times New Roman", size = 36L),
aspect.ratio = 1,
plot.margin = margin(10, 10, 10, 10))+
theme(legend.position = "none", #Xoá legend
plot.title = element_blank()) #Xoá title
# Hiển thị
print(F21)## `geom_smooth()` using formula = 'y ~ x'
#png("F21.png", width = 8, height = 8, units = "cm", res = 600)
#plot(F21, main = "My Base R Plot")
#dev.off()
#2222222222222222222222222222222222222222222222222222222222222222222222222222
library(ggplot2)
library(showtext)
showtext_auto()
model <- ket_qua_modelY2
# ========== TỰ ĐỘNG TÍNH GIỚI HẠN ĐÚNG NHƯ DESIGN-EXPERT ==========
n <- nrow(QHTN)
p <- length(coef(model))
df <- n - p - 1
alpha <- 0.05
critical_t <- qt(1 - alpha/(2*n), df) # Bonferroni critical value
lower_limit <- -critical_t
upper_limit <- +critical_t
cat(sprintf("Số quan sát n = %d, số tham số p = %d → df = %d\n", n, p, df))## Số quan sát n = 15, số tham số p = 8 → df = 6cat(sprintf("Giới hạn Bonferroni (α=0.05): ± %.4f\n", critical_t))## Giới hạn Bonferroni (α=0.05): ± 4.6979# Ví dụ: bạn sẽ thấy ±5.8, ±6.1, hoặc ±6.25 tùy model → khớp 100% với Design-Expert!
# =================================================================
esr <- rstudent(model)
run_order <- 1:length(esr)
df_run <- data.frame(
Run = run_order,
ESR = esr,
Y2 = QHTN$CiBD_Y2
#Y1 = QHTN$LucF_Y1
)
F22 <- ggplot(df_run, aes(x = Run, y = ESR)) +
geom_hline(yintercept = 0, color = "black", linewidth = 0.8) +
# Đường giới hạn TỰ ĐỘNG, CHÍNH XÁC NHƯ DESIGN-EXPERT
geom_hline(yintercept = lower_limit, color = "red", linetype = "solid", linewidth = 1) +
geom_hline(yintercept = upper_limit, color = "red", linetype = "solid", linewidth = 1) +
# Ghi chú giá trị giới hạn lên đồ thị (rất chuyên nghiệp)
annotate("text", x = 1, y = upper_limit,
label = sprintf("%.3f", upper_limit),
vjust = -0.8, hjust = 0, color = "red", size = 12, family = "Times New Roman") +
annotate("text", x = 1, y = lower_limit,
label = sprintf("%.3f", lower_limit),
vjust = 1.8, hjust = 0, color = "red", size = 12, family = "Times New Roman") +
geom_point(aes(fill = Y2), shape = 22, size = 2, color = "black", stroke = 0.5) +
geom_line(color = "black", linewidth = 0.5) +
scale_fill_gradientn(colours = viridis_col_simple,
name = expression(bold("Ci (Y2)"))) +
scale_x_continuous(name = "Run Number",
breaks = seq(1, max(run_order), by = 2), # chỉ hiện số lẻ: 1,3,5,...,15
minor_breaks = 1:max(run_order), # vẫn có gạch ngắn ở tất cả các số
expand = c(0.02, 0))+ # sát viền đẹp hơn
scale_y_continuous(name = "Externally Studentized Residuals",
limits = c(-6.2, 6.2),
#labels = function(x) sprintf("%.2f", x)
breaks = round(seq(-6, 6, 2), 1)) +
ggtitle("Residuals vs. Run Order\n(Externally Studentized Residuals with Bonferroni Limits)") +
theme_bw(base_family = "Times New Roman", base_size = 16) +
theme(
panel.grid.major.x = element_blank(), # Loại bỏ đường lưới dọc
panel.grid.minor.x = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.minor.y = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.major.y = element_blank(), # Loại bỏ đường lưới dọc
axis.ticks = element_line(linewidth = 0.5), #nét gạchight (tick marks) trên trục x và y
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5), # Thêm viền đen
axis.text = element_text(color = "black", size = 12), # Chữ trên trục màu đen
axis.title = element_text(color = "black", size = 12))+ # Tiêu đề trục màu đen
theme(plot.title = element_text(family = "Times New Roman Bold", size = 18L, hjust = 0.5),
axis.title.y = element_text(family = "Times New Roman", size = 36L, vjust = 2),
axis.title.x = element_text(family = "Times New Roman", size = 36L, , vjust = -2),
axis.text.y = element_text(family = "Times New Roman", size = 36L),
axis.text.x = element_text(family = "Times New Roman", size = 36L))+
theme(legend.title = element_text(family = "Times New Roman", size = 36L),
legend.text = element_text(family = "Times New Roman", size = 36L),
aspect.ratio = 1,
plot.margin = margin(10, 10, 10, 10))+
theme(legend.position = "none", #Xoá legend
plot.title = element_blank()) #Xoá title
print(F22)
#png("F22.png", width = 8, height = 8, units = "cm", res = 600)
#plot(F22, main = "My Base R Plot")
#dev.off()
#333333333333333333333333333333333333333333333333333333333333333333333333333
# Biểu đồ 3
# PREDICTED vs ACTUAL
# ===============================
library(ggplot2)
library(showtext)
showtext_auto()
model <- ket_qua_modelY2
#actual <- QHTN$LucF_Y1 # hoặc CiBD_Y2 nếu bạn đang làm Y2
actual <- QHTN$CiBD_Y2
predicted <- predict(model)
df_pred <- data.frame(
Actual = actual,
Predicted = predicted,
Y2 = actual # tô màu theo giá trị thực tế
)
# ---- Tự động xác định min–max và tạo breaks đẹp ----
range_val <- range(c(df_pred$Actual, df_pred$Predicted)) # min và max chung
min_val <- floor(range_val[1] / 0.1) * 0.1 # làm tròn xuống
max_val <- ceiling(range_val[2] / 0.1) * 0.1 # làm tròn lên
break_seq <- seq(min_val, max_val, by = 0.1) # bạn thay 5 → 10, 20… tùy thích
# ---- Biểu đồ hoàn hảo ----
F23 <- ggplot(df_pred, aes(x = Actual, y = Predicted)) +
# Đường hồi quy + vùng tin cậy 95%
geom_abline(intercept = 0, slope = 1, color = "black", linewidth = 0.5, linetype = "solid") +
# Các điểm tô màu giống contour(rsm)
geom_point(aes(fill = Y2), shape = 22, size = 2, color = "black", stroke = 0.5) +
# Palette terrain.colors giống hệt 2 đồ thị trước
scale_fill_gradientn(colours = viridis_col_simple,
name = expression(bold("Ci coefficient"))) +
# Trục x, y: tự động min–max, nhưng breaks tùy chỉnh
scale_x_continuous(name = "Actual",
limits = c(min_val, max_val),
breaks = break_seq) +
scale_y_continuous(name = "Predicted",
limits = c(min_val, max_val),
breaks = break_seq) +
ggtitle("Predicted vs. Actual Values") +
theme_bw(base_family = "Times New Roman", base_size = 12) +
theme(
panel.grid.major.x = element_blank(), # Loại bỏ đường lưới dọc
panel.grid.minor.x = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.minor.y = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.major.y = element_blank(), # Loại bỏ đường lưới dọc
axis.ticks = element_line(linewidth = 0.5), #nét gạchight (tick marks) trên trục x và y
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5), # Thêm viền đen
axis.text = element_text(color = "black", size = 12), # Chữ trên trục màu đen
axis.title = element_text(color = "black", size = 12))+ # Tiêu đề trục màu đen
theme(plot.title = element_text(family = "Times New Roman Bold", size = 12L, hjust = 0.5),
axis.title.y = element_text(family = "Times New Roman", size = 36L, vjust = 2),
axis.title.x = element_text(family = "Times New Roman", size = 36L, , vjust = -2),
axis.text.y = element_text(family = "Times New Roman", size = 36L),
axis.text.x = element_text(family = "Times New Roman", size = 36L),
aspect.ratio = 1,
plot.margin = margin(10, 10, 10, 10))+
theme(legend.position = "right", # hoặc "bottom" nếu muốn nằm dưới
legend.title = element_text(family = "Times New Roman", size = 36L, hjust = 0.5), # tiêu đề legend
legend.text = element_text(family = "Times New Roman", size = 36L), # số trên thanh màu
legend.key.height = unit(0.5, "cm"), # chiều cao thanh màu (tăng/giảm thoải mái)
legend.key.width = unit(0.5, "cm"), # độ dày thanh màu
#legend.background = element_rect(fill = "white", color = "black", linewidth = 0.5), # viền legend
legend.margin = margin(10, 10, 10, 10))+ # khoảng cách bên trong legend
theme(plot.title = element_blank()) #Xoá title
# Xem trước
print(F23)
#png("F23.png", width = 16, height = 8, units = "cm", res = 600)
#plot(F23, main = "My Base R Plot")
#dev.off()###3.2.2. Đồ thị One Factor
#4444444444444444444444444444444444444444444444444444444444444444444444444444
# Biểu đồ 4
# ===============================
# ONE FACTOR PLOT: Ảnh hưởng của e → Ci
# Trục Y tự động min–max + tùy chỉnh khoảng cách chia (ví dụ: mỗi 10)
# ===============================
library(ggplot2)
library(dplyr)
library(showtext)
showtext_auto()
model <- ket_qua_modelY2
# ---- 1. Tạo lưới dự báo cho e từ 15 → 35 mm (giữ x1 = x3 = 0) ----
R_real <- seq(15, 35, length.out = 200)
grid <- data.frame(
x1 = (R_real - 25)/10,
x2 = 0, # chuyển sang coded
x3 = 0
)
# Dự báo + 95% Confidence Interval
pred <- suppressWarnings(
predict(model, newdata = grid, interval = "confidence", level = 0.95)
)
grid_df <- data.frame(
R_real = R_real,
fit = pred[,"fit"],
lwr = pred[,"lwr"],
upr = pred[,"upr"]
)
# ---- 2. Lấy đúng 3 điểm trung tâm thực nghiệm (R = 5 mm) ----
center_points <- data %>%
filter(KCTam_e == 25 & GocLuon_R == 5) %>%
select(KCTam_e, CiBD_Y2)
# ---- 3. TỰ ĐỘNG xác định min–max trục Y + tùy chỉnh breaks ----
all_y_values <- c(grid_df$lwr, grid_df$upr, center_points$CiBD_Y2)
min_y <- floor(min(all_y_values) / 0.1) * 0.1 # làm tròn xuống đến bội số 10
max_y <- ceiling(max(all_y_values) / 0.1) * 0.1 # làm tròn lên đến bội số 10
# TỰ DO CHỈNH khoảng cách chia tại đây:
break_step <- 0.1 # ← Thay 10 → 20, 50, 5… tùy bạn
y_breaks <- seq((min_y)-0.1, max_y, by = break_step)
# ---- 4. Biểu đồ đẹp hơn cả Design-Expert ----
F24 <- ggplot() +
# Vùng 95% Confidence Interval
geom_ribbon(data = grid_df,
aes(x = R_real, ymin = lwr, ymax = upr),
fill = "gray", alpha = 0.3) +
# Đường dự báo trung bình
geom_line(data = grid_df,
aes(x = R_real, y = fit),
color = "black", linewidth = 0.5) +
# Chỉ 3 điểm đỏ ở tâm (R = 5)
geom_point(data = center_points,
aes(x = KCTam_e, y = CiBD_Y2),
color = "darkred", size = 1.2, shape = 19) +
# Trục X
scale_x_continuous(name = "Distance e (mm)",
breaks = seq(15, 35, by = 5)) +
# Trục Y: tự động min–max + bạn tùy chỉnh khoảng cách chia
scale_y_continuous(name = "Ci coefficient",
limits = c((min_y)-0.1, max_y),
breaks = y_breaks) + # ← Tùy chỉnh ở đây!
ggtitle("One Factor Plot: Effect of R on Ci\n(x1 = 0, x3 = 0)") +
# Theme đồng bộ tuyệt đối với tất cả biểu đồ trước
theme_bw(base_family = "Times New Roman", base_size = 12) +
theme(
panel.grid.major.x = element_blank(), # Loại bỏ đường lưới dọc
panel.grid.minor.x = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.minor.y = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.major.y = element_blank(), # Loại bỏ đường lưới dọc
axis.ticks = element_line(linewidth = 0.5), #nét gạchight (tick marks) trên trục x và y
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5), # Thêm viền đen
axis.text = element_text(color = "black", size = 36), # Chữ trên trục màu đen
axis.title = element_text(color = "black", size = 36))+ # Tiêu đề trục màu đen
theme(plot.title = element_text(family = "Times New Roman Bold", size = 18L, hjust = 0.5),
axis.title.y = element_text(family = "Times New Roman", size = 36L, vjust = 2),
axis.title.x = element_text(family = "Times New Roman", size = 36L, , vjust = -2),
axis.text.y = element_text(family = "Times New Roman", size = 36L),
axis.text.x = element_text(family = "Times New Roman", size = 36L))+
theme(legend.title = element_text(family = "Times New Roman", size = 36L),
legend.text = element_text(family = "Times New Roman", size = 36L),
aspect.ratio = 1,
plot.margin = margin(10, 5, 10, 5))+
theme(legend.position = "none", #Xoá legend
plot.title = element_blank()) #Xoá title
# Xem trước
print(F24)
#png("F24.png", width = 8, height = 8, units = "cm", res = 600)
#plot(F24, main = "My Base R Plot")
#dev.off()
#5555555555555555555555555555555555555555555555555555555555555555555555555555
# Biểu đồ 5
# ===============================
# ONE FACTOR PLOT: Ảnh hưởng của R → Ci
# Trục Y tự động min–max + tùy chỉnh khoảng cách chia (ví dụ: mỗi 10)
# ===============================
library(ggplot2)
library(dplyr)
library(showtext)
showtext_auto()
model <- ket_qua_modelY2
# ---- 1. Tạo lưới dự báo cho R từ 3 → 7 mm (giữ x1 = x3 = 0) ----
R_real <- seq(3, 7, length.out = 200)
grid <- data.frame(
x1 = 0,
x2 = (R_real - 5)/2, # chuyển sang coded
x3 = 0
)
# Dự báo + 95% Confidence Interval
pred <- suppressWarnings(
predict(model, newdata = grid, interval = "confidence", level = 0.95)
)
grid_df <- data.frame(
R_real = R_real,
fit = pred[,"fit"],
lwr = pred[,"lwr"],
upr = pred[,"upr"]
)
# ---- 2. Lấy đúng 3 điểm trung tâm thực nghiệm (R = 5 mm) ----
center_points <- data %>%
filter(GocLuon_R == 5 & KCTam_e == 25 ) %>%
select(GocLuon_R, CiBD_Y2)
# ---- 3. TỰ ĐỘNG xác định min–max trục Y + tùy chỉnh breaks ----
all_y_values <- c(grid_df$lwr, grid_df$upr, center_points$CiBD_Y2)
min_y <- floor(min(all_y_values) / 0.1) * 0.1 # làm tròn xuống đến bội số 10
max_y <- ceiling(max(all_y_values) / 0.1) * 0.1 # làm tròn lên đến bội số 10
# TỰ DO CHỈNH khoảng cách chia tại đây:
break_step <- 0.1 # ← Thay 10 → 20, 50, 5… tùy bạn
y_breaks <- seq(min_y, (max_y)+0.1, by = break_step)
# ---- 4. Biểu đồ đẹp hơn cả Design-Expert ----
F25 <- ggplot() +
# Vùng 95% Confidence Interval
geom_ribbon(data = grid_df,
aes(x = R_real, ymin = lwr, ymax = upr),
fill = "gray", alpha = 0.3) +
# Đường dự báo trung bình
geom_line(data = grid_df,
aes(x = R_real, y = fit),
color = "black", linewidth = 0.5) +
# Chỉ 3 điểm đỏ ở tâm (R = 5)
geom_point(data = center_points,
aes(x = GocLuon_R, y = CiBD_Y2),
color = "darkred", size = 1.2, shape = 19) +
# Trục X
scale_x_continuous(name = "Corner radius R (mm)",
breaks = seq(3, 7, by = 1)) +
# Trục Y: tự động min–max + bạn tùy chỉnh khoảng cách chia
scale_y_continuous(name = "Ci coefficient",
limits = c(min_y, (max_y)+0.1),
breaks = y_breaks) + # ← Tùy chỉnh ở đây!
ggtitle("One Factor Plot: Effect of R on Ci\n(x1 = 0, x3 = 0)") +
# Theme đồng bộ tuyệt đối với tất cả biểu đồ trước
theme_bw(base_family = "Times New Roman", base_size = 12) +
theme(
panel.grid.major.x = element_blank(), # Loại bỏ đường lưới dọc
panel.grid.minor.x = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.minor.y = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.major.y = element_blank(), # Loại bỏ đường lưới dọc
axis.ticks = element_line(linewidth = 0.5), #nét gạchight (tick marks) trên trục x và y
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5), # Thêm viền đen
axis.text = element_text(color = "black", size = 36), # Chữ trên trục màu đen
axis.title = element_text(color = "black", size = 36))+ # Tiêu đề trục màu đen
theme(plot.title = element_text(family = "Times New Roman Bold", size = 18L, hjust = 0.5),
axis.title.y = element_text(family = "Times New Roman", size = 36L, vjust = 2),
axis.title.x = element_text(family = "Times New Roman", size = 36L, , vjust = -2),
axis.text.y = element_text(family = "Times New Roman", size = 36L),
axis.text.x = element_text(family = "Times New Roman", size = 36L))+
theme(legend.title = element_text(family = "Times New Roman", size = 36L),
legend.text = element_text(family = "Times New Roman", size = 36L),
aspect.ratio = 1,
plot.margin = margin(10, 5, 10, 5))+
theme(legend.position = "none", #Xoá legend
plot.title = element_blank()) #Xoá title
# Xem trước
print(F25)
#png("F25.png", width = 8, height = 8, units = "cm", res = 600)
#plot(F25, main = "My Base R Plot")
#dev.off()
#6666666666666666666666666666666666666666666666666666666666666666666666666666
# Biểu đồ 6
# ===============================
# Biểu đồ One Factor: Ảnh hưởng của Hệ số ma sát μ → Ci
# Hoàn toàn giống format F4 của bạn (8x8 cm, font 72, không title/legend)
# ===============================
library(ggplot2)
library(dplyr)
library(showtext)
showtext_auto()
model <- ket_qua_modelY2 # hoặc ket_qua_model_opt
# ---- 1. Tạo lưới dự báo cho μ từ 0.10 → 0.20 (giữ x1=x2=0) ----
mu_real <- seq(0.10, 0.20, length.out = 200)
grid <- data.frame(
x1 = 0,
x2 = 0,
x3 = (mu_real - 0.15)/0.05 # chuyển sang coded: (μ - 0.15)/0.05
)
# Dự báo + 95% Confidence Interval
pred <- suppressWarnings(
predict(model, newdata = grid, interval = "confidence", level = 0.95)
)
grid_df <- data.frame(
mu_real = mu_real,
fit = pred[,"fit"],
lwr = pred[,"lwr"],
upr = pred[,"upr"]
)
# 3. LẤY ĐIỂM TRUNG TÂM – DÙNG round() ĐỂ TRÁNH LỖI LÀM TRÒN
center_points <- data %>%
mutate(mu_round = round(HesoMS_f, 3)) %>% # làm tròn 3 chữ số
filter(mu_round == 0.150 & KCTam_e == 25) %>% # chắc chắn bắt được 0.15
select(mu_real = HesoMS_f, CiBD_Y2) # đổi tên cột cho khớp
# Kiểm tra nhanh: phải ra đúng 3 dòng
print(center_points)## # A tibble: 3 × 2
## mu_real CiBD_Y2
## <dbl> <dbl>
## 1 0.15 0.458
## 2 0.15 0.483
## 3 0.15 0.509# ---- 3. TỰ ĐỘNG xác định min–max trục Y + tùy chỉnh breaks ----
all_y_values <- c(grid_df$lwr, grid_df$upr, center_points$CiBD_Y2)
min_y <- floor(min(all_y_values) / 0.1) * 0.1 # làm tròn xuống đến bội số 10
max_y <- ceiling(max(all_y_values) / 0.1) * 0.1 # làm tròn lên đến bội số 10
# TỰ DO CHỈNH khoảng cách chia tại đây:
break_step <- 0.1 # ← Thay 10 → 20, 50, 5…tùy bạn
y_breaks <- seq((min_y)-0.1, (max_y +0.1), by = break_step)
# ---- 4. Biểu đồ giống hệt F4 của bạn ----
F26 <- ggplot() +
geom_ribbon(data = grid_df,
aes(x = mu_real, ymin = lwr, ymax = upr),
fill = "gray", alpha = 0.3) +
geom_line(data = grid_df,
aes(x = mu_real, y = fit),
color = "black", linewidth = 0.5) +
geom_point(data = center_points,
aes(x = mu_real, y = CiBD_Y2),
color = "darkred", size = 1.2, shape = 19) +
scale_x_continuous(name = "Friction coefficient μ",
breaks = seq(0.10, 0.20, by = 0.02)) +
scale_y_continuous(name = "Ci coefficient",
limits = c((min_y)-0.1, (max_y +0.1)),
breaks = y_breaks) +
ggtitle("One Factor Plot: Effect of mu on Ci\n(x1 = 0, x3 = 0)") +
# Theme đồng bộ tuyệt đối với tất cả biểu đồ trước
theme_bw(base_family = "Times New Roman", base_size = 12) +
theme(
panel.grid.major.x = element_blank(), # Loại bỏ đường lưới dọc
panel.grid.minor.x = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.minor.y = element_blank(), # Loại bỏ đường lưới dọc nhỏ
panel.grid.major.y = element_blank(), # Loại bỏ đường lưới dọc
axis.ticks = element_line(linewidth = 0.5), #nét gạchight (tick marks) trên trục x và y
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5), # Thêm viền đen
axis.text = element_text(color = "black", size = 36), # Chữ trên trục màu đen
axis.title = element_text(color = "black", size = 36))+ # Tiêu đề trục màu đen
theme(plot.title = element_text(family = "Times New Roman Bold", size = 18L, hjust = 0.5),
axis.title.y = element_text(family = "Times New Roman", size = 36L, vjust = 2),
axis.title.x = element_text(family = "Times New Roman", size = 36L, , vjust = -2),
axis.text.y = element_text(family = "Times New Roman", size = 36L),
axis.text.x = element_text(family = "Times New Roman", size = 36L))+
theme(legend.title = element_text(family = "Times New Roman", size = 36L),
legend.text = element_text(family = "Times New Roman", size = 36L),
aspect.ratio = 1,
plot.margin = margin(10, 5, 10, 5))+
theme(legend.position = "none", #Xoá legend
plot.title = element_blank()) #Xoá title
# Xem trước
print(F26)
#png("F26.png", width = 8, height = 8, units = "cm", res = 600)
#plot(F26, main = "My Base R Plot")
#dev.off()###3.2.3. Đồ thị Contour Plot
#7777777777777777777777777777777777777777777777777777777777777777777777777777
## Biểu đồ contour ảnh hưởng e & R tới Ci
#png("F27v2.png", width = 8, height = 8, units = "cm", res = 600)
library(viridis)
# Tạo dải 120 màu (giống col_dx_120 của bạn)
viridis_col_120 <- viridis(120, option = "D") # Mặc định: viridis "D" (xanh tím → vàng) "turbo" là tím-đỏ
cividis_col_120 <- viridis(120, option = "E") # Xanh → vàng (colorblind-safe nhất)
my_col <- c("#0000FF", "#00FFFF", "#00FF00", "#FFFF00", "#FF8000", "#FF0000")
col_dx_120 <- colorRampPalette(my_col)(120)
viridis_col_full <- viridis(120, option = "turbo")
viridis_col_simple <- viridis_col_full[1:90] # Chọn từ vị trí 1 đến 90
par(family = "Times New Roman", ps = 36, cex = 1,
mar = c(4, 3, 1, 1), oma = c(0, 0, 0, 0), mgp = c(3, 0.6, 0), las = 1, tcl = -0.3)
# Tạo lưới đều
e_vals <- seq(15, 35, length.out = 300)
R_vals <- seq(3, 7, length.out = 300)
mu_vals <- seq(0.1, 0.2, length.out = 300)
# Mở rộng thêm 0.15 ở mỗi đầu → tạo đoạn thừa rõ rệt
xlim <- c(14.5, 35.5) # thay vì c(15,35)
ylim <- c(2.9, 7.1) # thay vì c(3,7)
grid <- expand.grid(x1 = (e_vals - 25)/10, x2 = (R_vals - 5)/2)
grid$x3 <- 0
grid$yhat <- predict(ket_qua_modelY2, grid)
zmat <- matrix(grid$yhat, 300, 300)
# Vẽ nền gradient đúng màu bạn muốn
image(e_vals, R_vals, zmat, col = viridis_col_simple, xlim = xlim, ylim = ylim,
xlab = "", ylab = "", axes = TRUE, useRaster = TRUE)
min_z <- min(zmat)
max_z <- max(zmat)
levels_05 <- seq(
floor(min_z/0.05)*0.05,
ceiling(max_z/0.05)*0.05,
by = 0.05)
# Vẽ đường đồng mức + nhãn số
contour(e_vals, R_vals, zmat, add = TRUE,
col = "black", lwd = 1, labcex = 1, drawlabels = TRUE, method = "flattest")
#levels = levels_05) # ★★★ ép mức đồng mức theo bước 0.05)
# Trục + nhãn
# Vẽ trục (chỉ vạch, không chữ)
#axis(1, tcl = -0.3, mgp = c(3, 0.6, 0))
#axis(2, tcl = -0.3, mgp = c(3, 0.8, 0), las = 1)
# THÊM NHÃN TRỤC ĐÚNG NHƯ BẠN MUỐN (đẹp, điều chỉnh được line, cex, font)
title(xlab = "Distance e (mm)",
line = 1.6, # chỉnh độ cao trục X ở đây
cex.lab = 1, # đúng 12 pt
family = "Times New Roman")
title(ylab = "Corner radius R (mm)",
line = 1.6, # chỉnh độ cao trục Y ở đây
cex.lab = 1, # đúng 12 pt
family = "Times New Roman")
# Tiêu đề slice (bạn tự đặt, không còn "Slice at x1 = 0")
mtext("Slice at Friction coefficient μ = 0.15",
family = "Times New Roman",
side = 1, line = 2.8, cex = 1, font = 1) # font = 2 → in đậm
#dev.off()###3.2.4. Đồ thị 3D Surface Plot
#===================================================================CHUẨN MẪU
#8888888888888888888888888888888888888888888888888888888888888888888888888888
## Biểu đồ 3D Surface ảnh hưởng e & R tới Ci
#png("F28v2.png", width = 8, height = 8, units = "cm", res = 600)
viridis_col_simple <- viridis(120, option = "D")
viridis_col_full <- viridis(120, option = "turbo")
viridis_col_simple <- viridis_col_full[1:90] # Chọn từ vị trí 1 đến 90
par(
family = "Times New Roman",
ps = 36,
cex = 1,
mar = c(2.6, 2.6, 0.1, 0.1), #(Dưới, Trái, Trên, Phải)
mgp = c(0.5, 0.02, 0.7) #(title, tick lable, tick mark)
)
# Tạo lưới đều
e_vals <- seq(15, 35, length.out = 300)
R_vals <- seq(3, 7, length.out = 300)
mu_vals <- seq(0.1, 0.2, length.out = 300)
grid <- expand.grid(x1 = (e_vals - 25)/10, x2 = (R_vals - 5)/2)
grid$x3 <- 0
grid$yhat <- predict(ket_qua_modelY2, grid)
zmat <- matrix(grid$yhat, 300, 300)
min_z <- min(zmat)
max_z <- max(zmat)
Color_Tren <- ((0.7-max_z)/(max_z - min_z))*90
Color_Duoi <-((min_z-0.2)/(max_z - min_z))*90
# Thêm các màu trống cả phía trước và phía sau
viridis_col_extended <- c(
rep(viridis_col_simple[1], Color_Duoi), # Thêm 36 màu giống màu đầu tiên (phía dưới)
viridis_col_simple, # Dải màu chính từ thấp đến cao
rep(viridis_col_simple[90], Color_Tren) # Thêm 36 màu giống màu cuối cùng (phía trên)
)
# Tạo các mức contour cách nhau 0.05
levels_05 <- seq(
floor(min_z/0.05)*0.05,
ceiling(max_z/0.05)*0.05,
by = 0.05)
# Vẽ đồ thị không có trục
res <- persp(
ket_qua_modelY2,
~ x1 + x2,
at = c(x3 = 0),
zlim = c(0.2, 0.7), # ← thêm khoảng trống !!!
axes = FALSE,
#sub = NULL, # ẩn Slice at
col = viridis_col_extended,
#contours = "colors",
contours = list(levels = levels_05, col = "colors"), # Mức bước 0.05, màu tự động
color.legend = TRUE,
expand = 1,
theta = -28,
phi = 20
)## Warning in persp.default(dat$x, dat$y, dat$z, zlim = dat$zlim, theta = theta, :
## "color.legend" is not a graphical parameter## Warning in persp.default(dat$x, dat$y, dat$z, xlab = dat$labs[1], ylab =
## dat$labs[2], : "color.legend" is not a graphical parameter## Warning in title(sub = dat$labs[5], ...): "axes" is not a graphical parameter## Warning in title(sub = dat$labs[5], ...): "color.legend" is not a graphical
## parameter## Warning in title(sub = dat$labs[5], ...): "expand" is not a graphical parameter# Tạo lưới và các giá trị tick
x1_seq_coded <- seq(min(QHTN$x1), max(QHTN$x1), length.out = 80)
x2_seq_coded <- seq(min(QHTN$x2), max(QHTN$x2), length.out = 80)
e_seq_coded <- seq(min(data$KCTam_e), max(data$KCTam_e), length.out = 80)
R_seq_coded <- seq(min(data$GocLuon_R), max(data$GocLuon_R), length.out = 80)
#f_seq_coded <- seq(min(data$HesoMS_f), max(data$HesoMS_f), length.out = 80)
grid <- expand.grid(x1 = x1_seq_coded, x2 = x2_seq_coded, x3 = 0)
grid_real <- expand.grid(x1 = e_seq_coded, x2 = R_seq_coded, x3 = 0.15)
grid$F_pred <- predict(ket_qua_modelY2, newdata = grid)
z_mat <- matrix(grid$F_pred, nrow = length(x1_seq_coded), ncol = length(x2_seq_coded))
x_ticks1 <- c("–","–","–","–","–")
x_ticks2 <- c("15","20","25","30","35")
y_ticks1 <- c("–","–","–","–","–")
y_ticks2 <- c("3","4","5","6","7")
#y_ticks1 <- c("–","–","–","–","–","–")
#y_ticks2 <- c("0.10","0.12","0.14","0.16","0.18","0.20")
z_ticks1 <- c("","–","–","–","–","–")
z_ticks2 <- c("","0.3","0.4","0.5","0.6","0.7")
par(xpd = NA)
#################################################### TRỤC X
###X11 Vẽ tick marks cho trục x2 với phép xoay toàn bộ trục
rotation_angle <- -73 # Góc xoay của toàn bộ trục x2 (độ)
angle_rad <- rotation_angle * pi / 180
offset1 <- 0.403
for(i in seq_along(x_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base11_x <- 0.247 #giảm là lên
base11_y <- -0.239 #giảm là sang trái
tick_offset <- (i/length(x_ticks1)) * offset1
x_before <- base11_x
y_before <- base11_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = x_ticks1[i],
cex = 1,
srt = rotation_angle)
}
###X12 Vẽ tick labels cho trục x2 với phép xoay toàn bộ trục
for(i in seq_along(x_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base12_x <- base11_x +0.03 #giảm là lên
base12_y <- base11_y #giảm là sang trái
tick_offset <- (i/length(x_ticks1)) * offset1
x_before <- base12_x
y_before <- base12_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = x_ticks2[i],
cex = 1,
srt = 0)
}
############ CÁCH 2: VẼ GẠCH TICK MARK BẰNG ĐIỂM
#for(i in seq_along(x_ticksmark)) {
# Tọa độ gốc của tick trước khi xoay
#tick_offset <- (i/length(x_ticksmark)) * 0.4
#x_before <- base_x
#y_before <- base_y + tick_offset
# Tính tọa độ sau khi xoay
#rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
#rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
#Vẽ tick mark tại vị trí đã xoay
#tmx21=rotated_x-0.062 #giảm là sang trái
#tmx22=rotated_y+0.01 #giảm là xuống
#segments(tmx21, tmx22,
#tmx21+0.008, tmx22-0.012,
#col = "black", lwd = 1.2)}
#################################################### TRỤC Y
###Y21 Vẽ tick marks cho trục x3 với phép xoay toàn bộ trục
rotation_angle <- 37 # Góc xoay của toàn bộ trục x3 (độ)
angle_rad <- rotation_angle * pi / 180
offset2 <- 0.245
for(i in seq_along(y_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base21_x <- -0.239 #giảm là sang trái
base21_y <- -0.22 #giảm là xuống
tick_offset <- (i/length(y_ticks1)) * offset2
x_before <- base21_x
y_before <- base21_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = y_ticks1[i],
cex = 1,
srt = rotation_angle)
}
###Y22 Vẽ tick labels cho trục x3 với phép xoay toàn bộ trục
for(i in seq_along(y_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base22_x <- base21_x -0.04 #giảm là sang trái
base22_y <- base21_y +0.01 #giảm là xuống
tick_offset <- (i/length(y_ticks1)) * offset2
x_before <- base22_x
y_before <- base22_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay CÁCH 2
#segments(rotated_x - 0.05, rotated_y,
#rotated_x + 0.05, rotated_y,
#col = "black", lwd = 1.2)
# Vẽ tick label tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = y_ticks2[i],
cex = 1,
srt = rotation_angle-37)
}
#################################################### TRỤC Z
###Z31 Vẽ tick marks cho trục z với phép xoay toàn bộ trục
rotation_angle <- 4.2 # Góc xoay của toàn bộ trục z (độ)
angle_rad <- rotation_angle * pi / 180
offset3 <- 0.3568
for(i in seq_along(z_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base31_x <- -0.217 #giảm là sang trái
base31_y <- -0.163 #giảm là xuống
tick_offset <- (i/length(z_ticks1)) * offset3
x_before <- base31_x
y_before <- base31_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = z_ticks1[i],
cex = 1,
srt = rotation_angle)
}
###Z32 Vẽ tick label cho trục z với phép xoay toàn bộ trục
for(i in seq_along(z_ticks1)) {
# Tọa độ gốc của tick trước khi xoay
base32_x <- base31_x -0.04 #giảm là sang trái
base32_y <- base31_y #giảm là xuống
tick_offset <- (i/length(z_ticks1)) * offset3
x_before <- base32_x
y_before <- base32_y + tick_offset
# Tính tọa độ sau khi xoay
rotated_x <- x_before * cos(angle_rad) - y_before * sin(angle_rad)
rotated_y <- x_before * sin(angle_rad) + y_before * cos(angle_rad)
# Vẽ tick mark tại vị trí đã xoay CÁCH 2
#segments(rotated_x - 0.05, rotated_y,
#rotated_x + 0.05, rotated_y,
#col = "black", lwd = 1.2)
# Vẽ tick label tại vị trí đã xoay
text(rotated_x, rotated_y,
labels = z_ticks2[i],
cex = 1,
srt = 0)
}
# Thêm nhãn cho các trục
par(xpd = NA, srt = -55)
text(-0.23, -0.24, "Corner radius R (mm)")
par(xpd = NA, srt = 20)
text(0.13, -0.3, "Distance e (mm)")
par(xpd = NA, srt = -85)
text(-0.32, 0.01, "Ci coefficient")
#dev.off()#4. Optimization (Tối ưu) ##4.1. Tính toán giá trị tối ưu
#==============================================================================
#********************************************************** TỐI ƯU HOÁ *****
#==============================================================================
# Multi-response optimization using desirability package
library(desirability)
# ----------------------------
# USER SETTINGS
# ----------------------------
goal_y1 <- "min" # LucF_Y1
goal_y2 <- "min" # CiBD_Y2
importance_y1 <- 3
importance_y2 <- 5
shape_y1 <- 1.0
shape_y2 <- 1.0
# Đảm bảo các giá trị thấp/cao không bị NA
y1_low <- min(QHTN$LucF_Y1); y1_target <- NA; y1_high <- max(QHTN$LucF_Y1)
y2_low <- 0.2; y2_target <- NA; y2_high <- max(QHTN$CiBD_Y2)
grid_n <- 11
top_start_n <- 6
# ----------------------------
if (!exists("ket_qua_model") || !exists("ket_qua_modelY2")) stop("Models not found")
x1_rng <- range(QHTN$x1); x2_rng <- range(QHTN$x2); x3_rng <- range(QHTN$x3)
e_min <- x1_rng[1]*10 + 25; e_max <- x1_rng[2]*10 + 25
R_min <- x2_rng[1]*2 + 5; R_max <- x2_rng[2]*2 + 5
mu_min <- x3_rng[1]*0.05 + 0.15; mu_max <- x3_rng[2]*0.05 + 0.15
cat(sprintf(" e: %.3f — %.3f\n R: %.3f — %.3f\n mu: %.4f — %.4f\n\n", e_min, e_max, R_min, R_max, mu_min, mu_max))## e: 15.000 — 35.000
## R: 3.000 — 7.000
## mu: 0.1000 — 0.2000# ----------------------------
# Define individual desirability functions
if (goal_y1 == "max") {
y1D <- dMax(y1_low, y1_high, scale = shape_y1)
} else if (goal_y1 == "min") {
y1D <- dMin(y1_low, y1_high, scale = shape_y1)
} else {
y1D <- dTarget(y1_low, y1_target, y1_high, scale = shape_y1)
}
if (goal_y2 == "max") {
y2D <- dMax(y2_low, y2_high, scale = shape_y2)
} else if (goal_y2 == "min") {
y2D <- dMin(y2_low, y2_high, scale = shape_y2)
} else {
y2D <- dTarget(y2_low, y2_target, y2_high, scale = shape_y2)
}
# LƯU Ý: Loại bỏ định nghĩa overallD do lỗi
cat("Using desirability bounds:\n")## Using desirability bounds:cat(sprintf(" Y1 (low, target, high) = %.4f, %s, %.4f\n", y1_low, ifelse(is.na(y1_target),"NA",format(y1_target, digits=6)), y1_high))## Y1 (low, target, high) = 107.1278, NA, 133.7514cat(sprintf(" Y2 (low, target, high) = %.4f, %s, %.4f\n\n", y2_low, ifelse(is.na(y2_target),"NA",format(y2_target, digits=6)), y2_high))## Y2 (low, target, high) = 0.2000, NA, 0.6149# ----------------------------
# Overall desirability function (TÍNH THỦ CÔNG)
overall_D_point <- function(par_real) {
e <- par_real[1]; R_real <- par_real[2]; mu_real <- par_real[3]
x1c <- (e - 25)/10; x2c <- (R_real - 5)/2; x3c <- (mu_real - 0.15)/0.05
newd <- data.frame(x1 = x1c, x2 = x2c, x3 = x3c)
y1 <- predict(ket_qua_model, newd)
y2 <- predict(ket_qua_modelY2, newd)
d1 <- predict(y1D, y1)
d2 <- predict(y2D, y2)
# Tính Desirability Tổng thể (D) theo công thức Geometric Mean
# D = (d1^w1 * d2^w2)^(1 / (w1 + w2))
D_val <- (d1^importance_y1 * d2^importance_y2)^(1 / (importance_y1 + importance_y2))
list(D = as.numeric(D_val), d1 = as.numeric(d1), d2 = as.numeric(d2), y1 = y1, y2 = y2,
coded = c(x1 = x1c, x2 = x2c, x3 = x3c), real = c(e = e, R = R_real, mu = mu_real))
}
# ----------------------------
# Coarse grid search (TÍNH D thủ công)
e_grid <- seq(e_min, e_max, length.out = grid_n)
R_grid <- seq(R_min, R_max, length.out = grid_n)
mu_grid <- seq(mu_min, mu_max, length.out = grid_n)
grid_coarse <- expand.grid(e = e_grid, R = R_grid, mu = mu_grid)
coded_coarse <- data.frame(x1 = (grid_coarse$e - 25)/10, x2 = (grid_coarse$R - 5)/2, x3 = (grid_coarse$mu - 0.15)/0.05)
grid_coarse$y1 <- predict(ket_qua_model, coded_coarse)
grid_coarse$y2 <- predict(ket_qua_modelY2, coded_coarse)
# Tính d1, d2 và D cho lưới thô
grid_coarse$d1 <- predict(y1D, grid_coarse$y1)
grid_coarse$d2 <- predict(y2D, grid_coarse$y2)
grid_coarse$D <- (grid_coarse$d1^importance_y1 * grid_coarse$d2^importance_y2)^(1 / (importance_y1 + importance_y2))
ord <- order(grid_coarse$D, decreasing = TRUE)
best_starts <- grid_coarse[ord[1:top_start_n], c("e","R","mu","y1","y2","d1","d2","D")]
cat("Top coarse candidates (real variables):\n")## Top coarse candidates (real variables):print(best_starts)## e R mu y1 y2 d1 d2 D
## 10 33 3 0.1 122.6922 0.2776648 0.4153889 0.8127966 0.6319158
## 11 35 3 0.1 122.5727 0.2800541 0.4198801 0.8070373 0.6316558
## 9 31 3 0.1 122.6392 0.2802323 0.4173812 0.8066079 0.6300338
## 8 29 3 0.1 122.4135 0.2877566 0.4258572 0.7884712 0.6258426
## 7 27 3 0.1 122.0153 0.3002378 0.4408167 0.7583865 0.6187688
## 6 25 3 0.1 121.4444 0.3176758 0.4622598 0.7163539 0.6078369# ----------------------------
# Optimization
obj_neg <- function(par_real) -overall_D_point(par_real)$D
best_list <- list()
for (i in 1:min(nrow(best_starts), top_start_n)) {
init <- as.numeric(best_starts[i, c("e","R","mu")])
lower <- c(e_min, R_min, mu_min); upper <- c(e_max, R_max, mu_max)
opt <- tryCatch(optim(init, obj_neg, method = "L-BFGS-B", lower = lower, upper = upper, control = list(maxit = 500)),
error = function(e) list(par = init, value = NA))
res_info <- overall_D_point(opt$par)
best_list[[i]] <- list(par = opt$par, info = res_info)
}
Dvals <- sapply(best_list, function(z) z$info$D)
best_idx <- which.max(Dvals)
best_final <- best_list[[best_idx]]
cat("\n=== BEST SOLUTION ===\n")##
## === BEST SOLUTION ===print(best_final$info)## $D
## [1] 0.6320734
##
## $d1
## [1] 0.4163485
##
## $d2
## [1] 0.811996
##
## $y1
## 1
## 122.6667
##
## $y2
## 1
## 0.2779969
##
## $coded
## x1 x2 x3
## 0.8768963 -1.0000000 -1.0000000
##
## $real
## e R mu
## 33.76896 3.00000 0.10000cat(sprintf("\nOverall D = %.6f\n", best_final$info$D))##
## Overall D = 0.632073cat("Real variables (e, R, mu):\n"); print(best_final$info$real)## Real variables (e, R, mu):## e R mu
## 33.76896 3.00000 0.10000cat("Coded variables (x1,x2,x3):\n"); print(best_final$info$coded)## Coded variables (x1,x2,x3):## x1 x2 x3
## 0.8768963 -1.0000000 -1.0000000cat("Predicted responses (y1,y2):\n"); print(round(c(best_final$info$y1, best_final$info$y2),6))## Predicted responses (y1,y2):## 1 1
## 122.666683 0.277997cat("Individual desirabilities (d1,d2):\n"); print(round(c(best_final$info$d1, best_final$info$d2),6))## Individual desirabilities (d1,d2):## [1] 0.416348 0.811996##4.2. Plot tối ưu
#===========================================================================
# VẼ BỘ BIỂU ĐỒ SOLUTION (Desirability Plots)
#===========================================================================
# Thêm thư viện cần thiết
library(ggplot2)
library(ggpubr) # Để dùng ggarrange
library(showtext)
# === 1. LẤY GIÁ TRỊ TỐI ƯU TRỰC TIẾP TỪ KẾT QUẢ TÍNH TOÁN ===
# Giả định: biến 'best_final' đã được tính toán ở phần code trước.
if (!exists("best_final")) {
stop("Lỗi: Biến 'best_final' (kết quả tối ưu) không tồn tại. Vui lòng chạy phần code tối ưu hóa trước.")
}
e_opt <- best_final$info$real['e']
R_opt <- best_final$info$real['R']
mu_opt <- best_final$info$real['mu']
F_opt <- best_final$info$y1
Ci_opt <- best_final$info$y2
D_opt <- best_final$info$D
# === 2. CÁC GIỚI HẠN ===
# Khoảng giới hạn của các biến/đáp ứng
e_min <- 15; e_max <- 35
R_min <- 3; R_max <- 7
mu_min <- 0.1; mu_max <- 0.2
F_min <- 107.128; F_max <- 133.751 # Khoảng giá trị thực nghiệm của Y1
Ci_min <- 0.2; Ci_max <- 0.614868 # Khoảng giá trị cho D-function của Y2
# TẠO HÀM VẼ TỪNG BIỂU ĐỒ (Dạng 'stepped' cho biến rời rạc hoặc 'sloped' cho đáp ứng)
#============================================================================
# ----------------------------
# USER SETTING: KÍCH THƯỚC CHỮ
# ---------------------------- *5.5 để export png 16x12 600dpi
SIZE_TEXT <- 4.3*2.5 # Điều chỉnh giá trị này để thay đổi cỡ chữ toàn bộ
# ----------------------------
# Thêm thư viện cần thiết
library(ggplot2)
library(ggpubr)
library(showtext)
# Đảm bảo đã load font Times New Roman
# === DỮ LIỆU TỐI ƯU (Giả định đã có) ===
if (!exists("best_final")) {
stop("Lỗi: Biến 'best_final' (kết quả tối ưu) không tồn tại. Vui lòng chạy phần code tối ưu hóa trước.")
}
e_opt <- best_final$info$real['e']; R_opt <- best_final$info$real['R']
mu_opt <- best_final$info$real['mu']; F_opt <- best_final$info$y1
Ci_opt <- best_final$info$y2; D_opt <- best_final$info$D
# Khoảng giới hạn của các biến/đáp ứng
e_min <- 15; e_max <- 35
R_min <- 3; R_max <- 7
mu_min <- 0.1; mu_max <- 0.2
F_min <- 107.128; F_max <- 133.751
Ci_min <- 0.2; Ci_max <- 0.614868
#111111111111111111111111111111111111111111111111111111111111111111111111111
# Hàm cho các yếu tố (A, B, C) - Dạng stepped
# ĐÃ SỬA: Thêm size_t vào tham số mặc định
plot_factor <- function(var_name, var_label, var_opt, var_min, var_max, size_t = SIZE_TEXT) {
Y_PADDING <- 1
var_range <- var_max - var_min
EXPANSION_FACTOR <- 0.1
var_min_line <- var_min - EXPANSION_FACTOR * var_range
var_max_line <- var_max + EXPANSION_FACTOR * var_range
Y_SCALE_COMPRESSION <- 0.8
Y_MAX_LIMIT <- 1.3
df_line_precise <- data.frame(
x = c(var_min_line, var_min, var_min, var_max, var_max, var_max_line),
y = c(0, 0, Y_SCALE_COMPRESSION, Y_SCALE_COMPRESSION, 0, 0)
)
ggplot() +
geom_hline(yintercept = -Y_PADDING, color = "gray40", linewidth = 0.5) +
geom_line(data = df_line_precise, aes(x = x, y = y),
color = "gray40", linewidth = 0.7) +
geom_point(aes(x = var_opt, y = Y_SCALE_COMPRESSION), color = "red", size = 2) +
annotate("text", x = var_min, y = 0, label = format(var_min, digits = 2),
hjust = 0, vjust = 1.5, family = "Times New Roman", size = size_t) +
annotate("text", x = var_max, y = 0, label = format(var_max, digits = 2),
hjust = 1, vjust = 1.5, family = "Times New Roman", size = size_t) +
annotate("text", x = mean(c(var_min, var_max)), y = -Y_PADDING / 2-0.1, # khoảng cách đến lề dưới khung Box
label = paste0(var_label, " = ", format(round(var_opt, 3))),
#format(resp_opt, digits = 6, nsmall = 4)
#digits = 3 là Tối thiểu $3$ chữ số có nghĩa
#nsmall = 4 là Tối thiểu $4$ chữ số sau dấu phẩy
hjust = 0.5, vjust = 0.5, family = "Times New Roman", size = size_t) +
scale_x_continuous(limits = c(var_min_line, var_max_line), expand = c(0, 0)) +
scale_y_continuous(limits = c(-Y_PADDING, Y_MAX_LIMIT), expand = c(0, 0)) +
labs(x = NULL, y = NULL) +
theme_minimal() +
theme(
plot.title = element_blank(), axis.title.x = element_blank(),
axis.title.y = element_blank(), axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.caption = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, linewidth=0.5),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
plot.margin = margin(5, 5, 5, 5, "pt")
)
}
#222222222222222222222222222222222222222222222222222222222222222222222222222
# Hàm cho các đáp ứng (Y1, Y2) - Dạng sloped
# ĐÃ SỬA: Thêm size_t vào tham số mặc định
plot_response <- function(resp_label, resp_opt, resp_min, resp_max, color_opt="blue", size_t = SIZE_TEXT) {
Y_PADDING <- 1 #Khoảng cách đường ngang dưới với lề dưới khung Box
resp_range <- resp_max - resp_min
EXPANSION_FACTOR <- 0.2
resp_min_line <- resp_min - EXPANSION_FACTOR * resp_range
resp_max_line <- resp_max + EXPANSION_FACTOR * resp_range
Y_MAX_LIMIT <- 1.6 #Khoảng cách đường ngang trên với lề trên khung Box
df_line_precise <- data.frame(
x = c(resp_min_line, resp_min, resp_max, resp_max_line),
y = c(1, 1, 0, 0)
)
x_rel_opt <- (resp_opt - resp_min) / resp_range
y_opt_pos <- 1 - x_rel_opt
ggplot() +
geom_hline(yintercept = -Y_PADDING, color = "gray40", linewidth = 0.5) +
geom_line(data = df_line_precise, aes(x = x, y = y),
color = "gray40", linewidth = 0.7) +
geom_point(aes(x = resp_opt, y = y_opt_pos), color = color_opt, size = 2) +
annotate("text", x = resp_min, y = 1, label = format(resp_min, digits = 3),
hjust = 1, vjust = -0.5, family = "Times New Roman", size = size_t) +
annotate("text", x = resp_max, y = 0, label = format(resp_max, digits = 3),
hjust = 0, vjust = -0.5, family = "Times New Roman", size = size_t) +
annotate("text", x = mean(c(resp_min, resp_max)), y = -Y_PADDING / 2 -0.1,
label = paste0(resp_label, " = ", format(resp_opt, digits = 6, nsmall = 4)),
hjust = 0.5, vjust = 0.5, family = "Times New Roman", size = size_t) +
scale_x_continuous(limits = c(resp_min_line, resp_max_line), expand = c(0, 0)) +
scale_y_continuous(limits = c(-Y_PADDING, Y_MAX_LIMIT), expand = c(0, 0)) +
labs(x = NULL, y = NULL) +
theme_minimal() +
theme(
plot.title = element_blank(), axis.title.x = element_blank(),
axis.title.y = element_blank(), axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.caption = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, linewidth=0.5),
panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
plot.margin = margin(5, 5, 5, 5, "pt")
)
}
# === GỌI HÀM VÀ GỘP BIỂU ĐỒ ===
plot_A <- plot_factor("A", "Distance e (mm)", e_opt, e_min, e_max)
plot_B <- plot_factor("B", "Corner radius R (mm)", R_opt, R_min, R_max)
plot_C <- plot_factor("C", "Friction coefficient μ", mu_opt, mu_min, mu_max)
plot_F <- plot_response("Force F (ton)", F_opt, F_min, F_max, color_opt="blue")
plot_Ci <- plot_response("Ci coefficient", Ci_opt, Ci_min, Ci_max, color_opt="blue")
# 6. TẠO BIỂU ĐỒ TRỐNG VÀ THÊM NHÃN DESIRABILITY VÀO BÊN TRONG BOX
plot_empty_with_desirability <- ggplot() +
theme_void() +
#theme(panel.border = element_rect(colour = "black", fill=NA, linewidth=0.5),
#plot.margin = margin(5, 5, 5, 5, "pt")) +
# Font và Size cho nhãn Desirability (Chính giữa Box)
annotate("text",
x = 0.5, y = 0.5,
label = paste0("Desirability = ", format(D_opt, digits=3),
"\nSolution 1 out of 57"),
hjust = 0.5, vjust = 0.5,
family = "Times New Roman", size = SIZE_TEXT,
lineheight = 1) #Khi export 16x12 600dpi thì chỉnh cái này
# TẠO KHUNG CHUNG VÀ GỘP BIỂU ĐỒ (ggarrange)
solution_plot <- ggpubr::ggarrange(
plot_A, plot_B,
plot_C, plot_F,
plot_Ci, plot_empty_with_desirability,
ncol = 2, nrow = 3,
align = "hv",
widths = c(1, 1),
heights = c(1, 1, 1),
common.legend = FALSE
)
# In hoặc lưu biểu đồ với độ rộng yêu cầu (ví dụ: 16cm)
print(solution_plot)
#png("solution2_plot.png", width = 16, height = 12, units = "cm", res = 600)
#plot(solution_plot, main = "My Base R Plot")
#dev.off()
# Nếu cần in ra code hoàn chỉnh cho môi trường R:
# cat(sprintf("Độ rộng mong muốn là 16. Sử dụng ggsave:\nggsave(\"solution_plot.png\", final_plot, width = 16, height = 15, units = \"cm\")\n"))