Air【省选模拟赛】

标签:计算几何

题目

这里写图片描述

分析

棱柱的体积正比于底面面积和三个支柱的高度和

把贡献的系数算出来,排序

贪心即可

code

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<queue>
#include<vector>
#include<set>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define dep(i,a,b) for(int i=a;i>=b;i--)
#define ll long long
#define mem(x,num) memset(x,num,sizeof x)
#define reg(x) for(int i=last[x];i;i=e[i].next)
#define lson (x<<1)
#define rson (x<<1|1)
#define mid ((l+r)>>1)
#define pi acos(-1)
#define mp make_pair
#define pb push_back
#define pa pair<int,int>
#define fi first
#define se second
#define lowbit(x) ((x)&(-x))
using namespace std;
inline ll read(){
    ll f=1,x=0;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
inline ll getgcd(ll x,ll y){return !x?y:getgcd(y%x,x);}
inline ll getlcm(ll x,ll y){return x/getgcd(x,y)*y;}
inline ll qpow(int x,int y,int p){
    ll re=1;
    while(y){
        if(y&1)re=(re*x)%p;
        y>>=1,x=(x*x)%p;
    }
    return re;
}
//**********head by yjjr**********
const int maxn=1e4+6;
struct point{double x,y;}p[maxn];
inline point operator + (point a,point b){return (point){a.x+b.x,a.y+b.y};}
inline point operator - (point a,point b){return (point){a.x-b.x,a.y-b.y};}
inline double operator * (point a,point b){return a.x*b.y-a.y*b.x;}
inline bool cmp(point a,point b){return atan2l(a.y,a.x)<atan2l(b.y,b.x);}
double val[maxn],h[maxn],ans;
int n,T;
int main()
{
    T=read();
    while(T--){
        n=read();mem(p,0);mem(h,0);mem(val,0);
        rep(i,1,n){
            int X=read(),Y=read();
            p[i].x=(double)X,p[i].y=(double)Y;
        }
        sort(p+1,p+1+n,cmp);
        p[++n]=p[1];p[0]=p[n-1];
        rep(i,1,n-1)val[0]+=p[i]*p[i+1]/6;
        rep(i,1,n-1)val[i]=(p[i-1]*p[i]+p[i]*p[i+1])/6;
        rep(i,0,n-1)h[i]=read();
        sort(val,val+n);sort(h,h+n);
        ans=0;
        rep(i,0,n-1)ans+=h[i]*val[i];
        printf("%.8lf\n",ans);
    }
    return 0;
}
%% 船舶航线优化主程序 clear; clc; close all; %% 参数设置 % 船舶参数(Suezmax型) ship.L_pp = 228.9; % 垂线间长(m) ship.B = 32.26; % 船宽(m) ship.T = 14.471; % 吃水(m) ship.Cb = 0.803; % 方形系数 ship.D = 8.5; % 螺旋桨直径(m) ship.sfoc = 104; % 燃油消耗率(g/kWh) ship.V = 12; % 设计航速(kn) ship.V_max = 16; % 最大航速(kn) ship.V_min = 10; % 最小航速(kn) ship.P_max = 32000; % 最大主机功率(kW) ship.rho_sea = 1026; % 海水密度(kg/m^3) ship.rho_air = 1.225; % 空气密度(kg/m^3) ship.A_T = 1918; % 水面上的最大横向面积(m^2) % 算法参数 population_size = 50; % 种群大小 max_generations = 100; % 最大迭代次数 pc = 0.8; % 交叉概率 pm = 0.1; % 变异概率 T0 = 100; % 初始温度 alpha = 0.95; % 退火系数 start_time = datetime(2025,2,12,18,0,0); % 航线参数 start_port = [35.2391, 119.7921]; % 日照港 end_port = [-28.5202, 114.3272]; % 杰拉尔顿港 max_segments = 30; % 最大航段数 lat_range = [-30, 55]; lon_range = [110, 160]; %% 气象数据加载与预处理 % 读取NetCDF文件 op_ncfile = 'C:\Users\dell\Downloads\41c65edbd5a653750e1d93b2a34fe6a8.nc'; wv_ncfile = 'C:\Users\dell\Downloads\abbf6737917f8802584d2b89ac9cf3ca.nc'; % 获取维度信息 op_lon = ncread(op_ncfile, 'longitude'); % 风场文件经度向量 op_lat = ncread(op_ncfile, 'latitude'); % 风场文件纬度向量 time = ncread(op_ncfile, 'valid_time'); % 时间(小时数,从1900-01-01起) wv_lon = ncread(wv_ncfile, 'longitude'); % 浪场文件经度向量 wv_lat = ncread(wv_ncfile, 'latitude'); % 浪场文件纬度向量 % 正确时间转换(基于NC文件元数据) time_ref = datetime(1970,1,1, 'TimeZone', 'UTC'); % 明确指定UTC时区 date_vec = time_ref + seconds(time); % 保持UTC时区 % 提取2025年1-2月数据 start_date = datetime(2025,2,12, 'TimeZone', 'UTC'); end_date = datetime(2025,2,27,23,59,59, 'TimeZone', 'UTC'); time_mask = (date_vec >= start_date) & (date_vec <= end_date); time_idx = find(time_mask); % 读取气象变量(优化维度处理) % 注意:NetCDF维度顺序通常为(longitude, latitude, time) u10 = permute(ncread(op_ncfile, 'u10', [1,1,1], [Inf, Inf, Inf]), [2,1,3]); % 调整为(lat,lon,time) v10 = permute(ncread(op_ncfile, 'v10', [1,1,1], [Inf, Inf, Inf]), [2,1,3]); swh = permute(ncread(wv_ncfile, 'swh', [1,1,1], [Inf, Inf, Inf]), [2,1,3]); cdww = permute(ncread(wv_ncfile, 'cdww', [1,1,1], [Inf, Inf, Inf]), [2,1,3]); % 创建纬度-经度网格 [LonGrid, LatGrid] = meshgrid(op_lon, op_lat); % 构建气象数据结构体 weather_data = struct(... 'op_lon', op_lon, ... % 经度向量 [nLon×1] 'op_lat', op_lat, ... % 纬度向量 [nLat×1] 'wv_lon', wv_lon, ... % 经度向量 [nLon×1] 'wv_lat', wv_lat, ... % 纬度向量 [nLat×1] 'LonGrid', LonGrid, ... % 经度网格 [nLat×nLon] 'LatGrid', LatGrid, ... % 纬度网格 [nLat×nLon] 'time', date_vec(time_idx), ... % 时间向量 [nTime×1] 'u10', u10(:,:,time_idx), ... % 风速U分量 [nLat×nLon×nTime] 'v10', v10(:,:,time_idx), ... % 风速V分量 [nLat×nLon×nTime] 'swh', swh(:,:,time_idx), ... % 有效波高 [nLat×nLon×nTime] 'cdww', cdww(:,:,time_idx)); % 波浪阻力系数 [nLat×nLon×nTime] %% 初始化种群 population = initialize_population(population_size, start_port, end_port, max_segments, lat_range, lon_range); %% 混合遗传算法主循环 best_fitness = inf; best_route = []; best_fuel = Inf; best_time = Inf; fitness_history = zeros(max_generations, 1); temperature = T0; for gen = 1:max_generations % 计算适应度 fitness = zeros(population_size, 1); fuel_consumptions = zeros(population_size, 1); total_times = zeros(population_size, 1); for i = 1:population_size [fitness(i), fuel_consumptions(i), total_times(i)] = calculate_fitness(population{i}, ship, weather_data, start_time); end % 记录最佳个体 [min_fitness, idx] = min(fitness); current_fuel = fuel_consumptions(idx); current_time = total_times(idx); if min_fitness < best_fitness best_fitness = min_fitness; best_route = population{idx}; best_fuel = current_fuel; best_time = current_time; end fitness_history(gen) = best_fitness; % 择操作 selected = tournament_selection(population, fitness, population_size/2); % 交叉操作 offspring = hybrid_crossover(selected, pc, start_port, end_port); % 变异操作 offspring = mutation(offspring, pm, lat_range, lon_range); % 合并种群 combined_population = [population; offspring]; offspring_fitness = zeros(length(offspring), 1); for i = 1:length(offspring) offspring_fitness(i) = calculate_fitness(offspring{i}, ship, weather_data, start_time); end combined_fitness = [fitness; offspring_fitness]; % 模拟退火择 population = simulated_annealing_selection(combined_population, combined_fitness, population_size, temperature); % 更新温度 temperature = alpha * temperature; % 显示进度 fprintf('Generation %d: Fitness=%.2f, Fuel=%.2fkg, Time=%.2fh (Current) | Best: Fitness=%.2f, Fuel=%.2fkg, Time=%.2fh\n',... gen, min_fitness, current_fuel, current_time, best_fitness, best_fuel, best_time); % 更新温度 temperature = alpha * temperature; end %% 结果可视化 plot_route(best_route, start_port, end_port, weather_data); figure; plot(1:max_generations, fitness_history); xlabel('Generation'); ylabel('Best Fitness'); title('Convergence Curve'); %% 初始化种群函数 function population = initialize_population(pop_size, start_port, end_port, max_segments, lat_range, lon_range) population = cell(pop_size, 1); % 初始化为列向量 for i = 1:pop_size num_segments = randi([2, max_segments]); route = zeros(num_segments+1, 2); route(1,:) = start_port; route(end,:) = end_port; % 生成中间点 for j = 2:num_segments while true new_lat = lat_range(1) + (lat_range(2)-lat_range(1))*rand(); new_lon = lon_range(1) + (lon_range(2)-lon_range(1))*rand(); if ~isequal([new_lat, new_lon], start_port) && ~isequal([new_lat, new_lon], end_port) route(j,:) = [new_lat, new_lon]; break; end end end population{i} = route; % 按经度排序中间点,确保东向航行 route(2:end-1,:) = sortrows(route(2:end-1,:), 2); end end %% 适应度计算函数 function [fitness, total_fuel, total_time] = calculate_fitness(route, ship, weather_data, start_time) % 检查路径有效性 if isempty(route) || size(route,1) < 2 || any(isnan(route(:))) fitness = Inf; return; end total_distance = 0; total_time = 0; total_fuel = 0; % 计算每个航段的参数 for i = 1:size(route,1)-1 % 计算航段距离 (大圆距离) [dist, ~] = lldistkm(route(i,:), route(i+1,:)); if isnan(dist) % 如果距离计算失败 fitness = Inf; return; end total_distance = total_distance + dist; % 获取航段中点位置和时间的气象数据 mid_point = [(route(i,1)+route(i+1,1))/2, (route(i,2)+route(i+1,2))/2]; segment_time = start_time + seconds(total_time*3600); segment_time.TimeZone = 'UTC'; % 强制设置时区 % 双线性插值获取气象数据 [u10, v10, swh, cdww] = get_weather_data(mid_point, segment_time, weather_data); % 计算航向角 (0-360度,正北为0) heading = calculate_heading(route(i,:), route(i+1,:)); % 计算相对风向和波向 wind_dir = atan2d(v10, u10); % 风向(度) wind_speed = sqrt(u10^2 + v10^2); rel_wind_angle = wrapTo360(wind_dir - heading); % 计算船舶实际速度 (考虑风浪影响) [V_actual, fuel_consumption] = calculate_ship_performance(ship, wind_speed, rel_wind_angle, swh, cdww); % 计算航段时间和燃油消耗 segment_time_hours = dist / (V_actual * 1.852); % 海里转公里 total_time = total_time + segment_time_hours; total_fuel = total_fuel + fuel_consumption * segment_time_hours; % 添加航速有效性检查 if V_actual <= 0 || isnan(V_actual) fitness = Inf; return; end % 添加时间计算保护 segment_time_hours = max(0, dist / (V_actual * 1.852)); end % 适应度函数 (加权和) fuel_weight = 0.6; % 燃油消耗权重 fitness = fuel_weight * total_fuel + (1-fuel_weight) * total_time; fprintf('Segment %d: V=%.2fkn, Fuel=%.2fkg/h, Time=%.2fh\n',... i, V_actual, fuel_consumption, segment_time); end %% 气象数据获取函数 function [u10, v10, swh, cdww] = get_weather_data(point, time, weather_data) % 双线性插值获取气象数据 % 风场数据 (0.25°分辨率) [u10, v10] = bilinear_interp(point(1), point(2), time, ... weather_data.op_lat, weather_data.op_lon, weather_data.time, ... weather_data.u10, weather_data.v10); % 浪场数据 (0.5°分辨率) [swh, cdww] = bilinear_interp(point(1), point(2), time, ... weather_data.wv_lat, weather_data.wv_lon, weather_data.time, ... weather_data.swh, weather_data.cdww); end %% 双线性插值法函数 function [varargout] = bilinear_interp(lat, lon, time, lat_grid, lon_grid, time_grid, varargin) % 通用双线性插值函数 varargout = cell(1, nargin-6); % 找到最近的4个空间点和2个时间点 [lat_idx, lon_idx, time_idx, weights] = find_nearest_points(lat, lon, time, lat_grid, lon_grid, time_grid); for i = 1:length(varargin) var = varargin{i}; % 空间插值 spatial_interp = (1-weights(1))*(1-weights(2)) * var(lat_idx(1), lon_idx(1), time_idx(1)) + ... (1-weights(1))*weights(2) * var(lat_idx(1), lon_idx(2), time_idx(1)) + ... weights(1)*(1-weights(2)) * var(lat_idx(2), lon_idx(1), time_idx(1)) + ... weights(1)*weights(2) * var(lat_idx(2), lon_idx(2), time_idx(1)); % 时间插值 if length(time_idx) > 1 spatial_interp2 = (1-weights(1))*(1-weights(2)) * var(lat_idx(1), lon_idx(1), time_idx(2)) + ... (1-weights(1))*weights(2) * var(lat_idx(1), lon_idx(2), time_idx(2)) + ... weights(1)*(1-weights(2)) * var(lat_idx(2), lon_idx(1), time_idx(2)) + ... weights(1)*weights(2) * var(lat_idx(2), lon_idx(2), time_idx(2)); varargout{i} = (1-weights(3)) * spatial_interp + weights(3) * spatial_interp2; else varargout{i} = spatial_interp; end end end %% 最近的网格点寻找函数 function [lat_idx, lon_idx, time_idx, weights] = find_nearest_points(lat, lon, time, lat_grid, lon_grid, time_grid) % 找到最近的网格点和插值权重 % 纬度 [~, lat_idx(1)] = min(abs(lat_grid - lat)); if lat_grid(lat_idx(1)) < lat lat_idx(2) = min(lat_idx(1)+1, length(lat_grid)); else lat_idx(2) = max(lat_idx(1)-1, 1); end weights(1) = (lat - lat_grid(lat_idx(1))) / (lat_grid(lat_idx(2)) - lat_grid(lat_idx(1))); % 经度 [~, lon_idx(1)] = min(abs(lon_grid - lon)); if lon_grid(lon_idx(1)) < lon lon_idx(2) = min(lon_idx(1)+1, length(lon_grid)); else lon_idx(2) = max(lon_idx(1)-1, 1); end weights(2) = (lon - lon_grid(lon_idx(1))) / (lon_grid(lon_idx(2)) - lon_grid(lon_idx(1))); % 时间 time_diff = abs(time_grid - time); [~, time_idx(1)] = min(time_diff); if time_grid(time_idx(1)) < time time_idx(2) = min(time_idx(1)+1, length(time_grid)); else time_idx(2) = max(time_idx(1)-1, 1); end if length(time_idx) > 1 weights(3) = seconds(time - time_grid(time_idx(1))) / seconds(time_grid(time_idx(2)) - time_grid(time_idx(1))); end end %% 角度计算函数 function heading = calculate_heading(point1, point2) % 计算两点之间的航向角 (0-360度,正北为0) lat1 = deg2rad(point1(1)); lon1 = deg2rad(point1(2)); lat2 = deg2rad(point2(1)); lon2 = deg2rad(point2(2)); dLon = lon2 - lon1; y = sin(dLon) * cos(lat2); x = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dLon); heading = wrapTo360(rad2deg(atan2(y, x))); end %% 船舶油耗及速度计算函数 function [V_actual, fuel_consumption] = calculate_ship_performance(ship, wind_speed, rel_wind_angle, swh, cdww) % 计算考虑风浪影响的船舶实际速度和燃油消耗 % 计算风阻力 R_wind = 0.5 * ship.rho_air * wind_speed^2 * ship.A_T * abs(cosd(rel_wind_angle)); % 计算波浪阻力 (简化模型) R_wave = 0.5 * ship.rho_sea * (swh/2)^2 * ship.B * ship.T * (1 + 4 * (ship.T/ship.L_pp)^2); % 计算静水阻力 R_water = 0.5 * ship.rho_sea * ship.L_pp * ship.B * cdww; V_knots = ship.V; % 初始使用设计航速 V_ms = V_knots * 0.5144; % 节转m/s % 总静水阻力 R_total = R_wind + R_wave + R_water; % 考虑风浪影响的速度损失 (经验公式) speed_loss = 0.05 * (swh/4) + 0.02 * wind_speed/10; V_actual = max(ship.V_min, ship.V * (1 - speed_loss)); % 计算燃油消耗 (简化模型) V_actual_ms = V_actual * 0.5144; % 节转m/s power_required = R_total * V_actual_ms / 0.6; % 推进效率0.6 fuel_consumption = power_required * ship.sfoc / 1e6; % kg/h end %% 择算子 function selected = tournament_selection(population, fitness, num_selected) % 锦标赛择 selected = cell(num_selected, 1); for i = 1:num_selected % 随机择2个个体进行竞争 candidates = randperm(length(population), 2); if fitness(candidates(1)) < fitness(candidates(2)) selected{i} = population{candidates(1)}; else selected{i} = population{candidates(2)}; end end end %% 交叉算子 function offspring = hybrid_crossover(parents, pc, start_port, end_port) offspring = cell(length(parents), 1); % 确保子代数量与父代相同 i = 1; while i <= length(parents) if i+1 > length(parents) % 处理奇数情况:最后一个父代直接复制 offspring{i} = parents{i}; break; end if rand() < pc parent1 = parents{i}; parent2 = parents{i+1}; % 检查父代路径有效性 if size(parent1,1) < 3 || size(parent2,1) < 3 offspring{i} = parent1; offspring{i+1} = parent2; i = i + 2; continue; end min_len = min(size(parent1,1), size(parent2,1)); if min_len <= 2 offspring{i} = parent1; offspring{i+1} = parent2; i = i + 2; continue; end % 随机交叉点 crossover_point = randi([2, min_len-1]); % 生成子代 child1 = [parent1(1:crossover_point,:); parent2(crossover_point+1:end,:)]; child2 = [parent2(1:crossover_point,:); parent1(crossover_point+1:end,:)]; % 确保子代有效性 child1(1,:) = start_port; child1(end,:) = end_port; child2(1,:) = start_port; child2(end,:) = end_port; offspring{i} = child1; offspring{i+1} = child2; else % 不交叉时直接复制父代 offspring{i} = parents{i}; offspring{i+1} = parents{i+1}; end i = i + 2; % 移动到下一对父代 end end %% 变异算子 function mutated_pop = mutation(population, pm, lat_range, lon_range) mutated_pop = cell(size(population)); for i = 1:length(population) if rand() < pm route = population{i}; if size(route,1) <= 2 mutated_pop{i} = route; continue; end idx = randi([2, size(route,1)-1]); % 随机择中间点 new_lat = lat_range(1) + (lat_range(2)-lat_range(1))*rand(); new_lon = lon_range(1) + (lon_range(2)-lon_range(1))*rand(); route(idx,:) = [new_lat, new_lon]; mutated_pop{i} = route; else mutated_pop{i} = population{i}; end end end %% 模拟退火 function new_population = simulated_annealing_selection(population, fitness, pop_size, temperature) % 模拟退火择 [~, sorted_idx] = sort(fitness); new_population = cell(pop_size, 1); % 保留最优个体 new_population{1} = population{sorted_idx(1)}; % 对其他个体进行退火择 for i = 2:pop_size % 随机择两个个体 candidates = randperm(length(population), 2); delta = fitness(candidates(2)) - fitness(candidates(1)); % 接受概率 if delta < 0 new_population{i} = population{candidates(2)}; else p = exp(-delta / temperature); if rand() < p new_population{i} = population{candidates(2)}; else new_population{i} = population{candidates(1)}; end end end end %% 绘图函数 function plot_route(route, start_port, end_port, weather_data) % 绘制优化后的航线 figure; % 绘制地图背景 worldmap([-40 60], [100 170]); load coastlines; plotm(coastlat, coastlon); % 绘制航线 plotm(route(:,1), route(:,2), 'r-o', 'LineWidth', 2); plotm(start_port(1), start_port(2), 'go', 'MarkerSize', 10, 'MarkerFaceColor', 'g'); plotm(end_port(1), end_port(2), 'bo', 'MarkerSize', 10, 'MarkerFaceColor', 'b'); % 添加标题和图例 title('Optimized Ship Route from Dalian to Gladstone'); legend('Route', 'Dalian', 'Gladstone', 'Location', 'best'); % 添加气象数据示例 (可) hold on; [LON, LAT] = meshgrid(weather_data.op_lon, weather_data.op_lat); quiverm(LAT(1:5:end,1:5:end), LON(1:5:end,1:5:end), ... weather_data.u10(1:5:end,1:5:end,1), weather_data.v10(1:5:end,1:5:end,1), 'k'); hold off; end %% 距离及角度计算函数 function [km, arc] = lldistkm(latlon1, latlon2) % 计算两点之间的大圆距离 lat1 = latlon1(1); lon1 = latlon1(2); lat2 = latlon2(1); lon2 = latlon2(2); R = 6371; % 地球半径(km) lat1 = deg2rad(lat1); lon1 = deg2rad(lon1); lat2 = deg2rad(lat2); lon2 = deg2rad(lon2); delta_lat = lat2 - lat1; delta_lon = lon2 - lon1; a = sin(delta_lat/2)^2 + cos(lat1) * cos(lat2) * sin(delta_lon/2)^2; c = 2 * atan2(sqrt(a), sqrt(1-a)); km = R * c; arc = rad2deg(c); end 那你把这个代码修改一下发给我
最新发布
05-28
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值